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Abstract 


Title  of  Thesis:  Smart  Structures  and  Wavelet  Based 
System  Identification 
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In  this  thesis  we  study  the  problem  of  modeling  and  control  of  vibrations 
of  a  flexible  structure.  Different  approaches  have  been  considered  for  fitting  a 
finite  dimensional  model  to  the  infinite  dimensional  system.  In  particular  we 
used  a  discretization  method  to  model  the  cantilever  beam  and  using  this  model 
and  taking  into  account  the  saturation  of  input  signal,  we  have  designed  and 
analyzed  a  nonlinear  controller  .  An  observer  based  version  of  the  controller  has 
also  been  implemented  and  shown  to  be  stable. 

A  simulation  program  is  written  for  solving  the  governing  partial  differential 
equation  for  a  cantilever  beam  and  the  effect  of  passive  damping  due  to  viscous 
and  internal  damping  has  been  considered  to  get  a  more  realistic  simulation. 

Finally,  a  general  algorithm  for  nonparametric  system  identification  is  devel¬ 
oped  and  implemented  and  is  used  to  obtain  an  approximation  to  the  transfer 


function  matrix  of  a  smart  structure.  The  results  are  compared  with 
tional  methods. 
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Chapter  1 


Introduction 


Motivated  by  the  rapid  developments  of  problems  of  design  of  large  flexible 
space  structures,  vibration  suppression  by  active  modal  control  has  drawn  much 
attention  in  recent  years  and  has  been  the  subject  of  much  research. 

Generally,  large  spacecraft /space  structures  have  very  poor  passive  damping 
due  to  the  materials  used  and  also  the  fact  that  there  might  be  no  air  to  induce 
viscous  damping.  Another  characteristic  of  these  space  structures  is  that  the 
natural  frequency  of  the  structure  (at  least  the  lowest)  may  fall  into  the  band¬ 
width  of  the  attitude  control  system  resulting  in  large  amplitude  vibrations  and 
fatigue  in  the  whole  structure. 

Mechanical  systems  of  this  type  are  generally  governed  by  a  system  of  partial 
differential  equations;  thus,  evolution  equations  of  them  may  be  modeled  in 
infinite  dimensional  state  space. 

Despite  the  recent  advances  in  the  theory  of  distributed-parameter  control, 
(see  [19]  and  [14]  for  example),  it  is  still  a  growing  field.  Also,  usually 
distributed-parameter  control  requires  some  sort  of  distributed  sensor  which 
may  be  costly  and  difficult  to  implement.  As  a  result,  there  is  always  a  tendency 
toward  the  approximation  of  the  infinite  dimensional  system  with  a  finite  dimen- 
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sional  one  to  take  advantage  of  the  existing  well-developed  lumped-parameter 
control  theory.  Of  course,  this  approach  (truncation/ approximation)  may  entail 
some  performance  trade-off. 

There  has  been  a  significant  effort  in  the  field  of  distributed  parameter  control 
with  specific  application  in  active  vibration  damping  of  flexible  structures,  e.g. 
Balas  [11]  devised  a  control  scheme  for  a  certain  class  of  flexible  structures  and 
Curtain  and  Glover  [16]  addresses  the  problem  of  designing  a  finite  dimensional 
controller  for  the  control  of  a  distributed-parameter  system. 

The  major  goal  in  the  first  two  chapters  (chapters  two  and  three)  of  this 
thesis  is  to  exploit  some  kind  of  finite  dimensional  approximation  method  to 
obtain  a  model  (an  approximation)  for  a  cantilever  beam.  In  chapter  two,  a 
Galerkin  approximation  method  is  used  to  model  the  beam  with  several  rigid 
rods  that  are  compliantly  linked  together,  and  equations  of  motion  are  derived 
based  on  this  model.  The  rest  of  the  chapter  is  devoted  to  obtaining  estimates  for 
various  parameters  that  are  involved  in  the  model.  Once  the  system  of  ordinary 
differential  equations  has  been  derived,  there  is  another  step  of  simplification 
that  is  needed,  namely  linearization. 

The  detailed  models  of  piezoceramic  actuated  beam  can  be  quite  involved 
and  complicated  (see  [4]  for  example).  Usually  several  aspects  such  as  stiff¬ 
ness  of  the  bonding  material  or  the  way  in  which  the  piezo  ceramic  is  mounted 
(e.g.  embedded  or  bonded)  comes  into  the  picture.  For  example,  in  a  real  world 
application  like  active  vibration  damping  in  a  helicopter  rotor  [17],  the  heli¬ 
copter  blade  cannot  be  precisely  modeled  by  an  Euler-Bernoulli  beam  and  more 
complicated  models  should  be  used  that  in  turn  need  the  knowledge  of  many 
physical  parameters  that  are  not  so  easily  obtainable. 

In  chapter  three  we  developed  a  general  non-parametric  system  identification 
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tool  that  uses  the  idea  of  rational  wavelet  decomposition.  The  matching  pursuit 
algorithm  of  Mall  at  et.  ah  [18]  has  been  applied  to  the  problem  of  approximation 
of  a  given  frequency  response  with  a  real  rational  transfer  function.  This  method 
is  then  used  to  obtain  a  finite  dimensional  model  of  the  experimental  set  up  that 
consists  of  a  flexible  beam  and  piezoceramic  actuators. 

Chapter  four  addresses  the  problem  of  control.  First,  a  classical  LQR  type 
controller  is  employed  to  control  the  vibration.  This  is  later  used  as  a  basis 
for  comparison  with  other  methods.  In  a  real  control  set-up,  however,  there 
is  a  limit  on  the  amplitude  of  the  voltage  that  can  be  applied  safely  to  the 
actuating  piezoceramics  and  this  saturation  limit  on  control  signals  may  lead 
to  instability  or  poor  performance  of  the  system  if  one  does  not  take  voltage 
limit  into  account  in  the  design  of  the  controller.  Stability  of  systems  with 
bounded  controls  is  a  topic  studied  in  depth  by  Sussman  et  al.  [21],  Gutman 
and  Hagander  [13]  and  others.  Gutman  and  Hagander  proposed  a  control 
algorithm  to  accommodate  bounded  controls  in  [13].  In  chapter  four,  a  controller 
of  this  type  has  been  designed  and  tested  in  a  simulation  program.  This  type  of 
controller,  as  presented  in  [13]  needs  state  variable  information  to  implement  the 
stabilization,  a  requirement  that  may  not  be  feasible  in  most  practical  situations. 

In  chapter  five,  we  provide  a  modification  to  the  Gutman-Hagander  algo¬ 
rithm  to  allow  for  the  usage  of  estimated  state  variables  obtained  by  a  classical 
state  observer.  The  proof  of  stability  is  also  given  in  chapter  five.  Although  the 
main  goal  in  this  type  of  controller  is  stability,  an  improvement  in  settling  time 
of  the  response  is  observed  in  comparison  to  the  LQR  method. 

Chapter  six  offers  an  algorithm  (implicit  finite  difference)  for  finding  the  nu¬ 
merical  solution  to  the  PDE  describing  the  behavior  of  a  cantilever  beam  with 
piezo  actuator  material  mounted  on  it  and  piezoceramics  used  as  sensors.  The 
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effect  of  air  drag  and  internal  damping  is  also  included  in  the  model.  This  can 
serve  as  a  simulation  tool  to  test  the  performance  of  control  strategies.  With 
the  application  of  spectral  analysis  method  and  the  simulation  program  at  hand, 
the  empirical  frequency  response  of  the  cantilever  beam  with  piezoceramic  ac¬ 
tuators  and  sensors  is  obtained  and  is  later  used  for  the  purpose  of  system 
identification.  The  final  section  of  chapter  six  has  been  devoted  to  the  deriva¬ 
tion  of  an  analytical  expression  for  the  transfer  function  of  a  cantilever  beam 
with  bonded  piezoceramics.  The  empirical  frequency  response  obtained  through 
spectral  analysis  is  shown  to  corroborate  the  one  obtained  from  the  closed  form 
representation  of  the  transfer  function. 

Chapter  seven  contains  results,  some  concluding  remarks,  and  suggestions 
for  future  work  in  this  area. 
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Chapter  2 


Galerkin  Model  of  a  Flexible 
Cantilever  Beam 


First,  equations  of  motion  of  an  approximated  finite  dimensional  model  for  a 
cantilever  beam  with  embedded  piezo  actuators  are  obtained.  The  model  con¬ 
sists  of  N  rigid  bodies  that  are  linked  together  via  torsional  springs.  We  use 
piezoceramics  as  actuators  which  apply  point  moments  at  the  joints.  These 
actuators  may  have  non-uniform  spatial  distribution.  Based  on  the  physical  pa¬ 
rameters  of  the  beam,  some  expressions  for  the  stiffness  coefficients  and  moment 
generators  that  model  the  effect  of  piezoceramics  has  been  obtained.  Finally, 
simulations  are  also  performed  to  verify  the  correctness  of  the  model. 


2.1  Derivation  of  Equations  of  Motion 

A  flexible  beam  is  an  infinite  dimensional  dynamical  system.  In  order  to  use 
the  classical  control  methods  for  damping  the  vibration  of  a  flexible  cantilever 
beam,  we  first  need  to  approximate  it  with  a  finite  dimensional  model. 

Consider  the  following  model  in  which  the  beam  is  approximated  with  N 
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rigid  links,  connected  to  each  other  by  torsional  springs.  Figure  2.1  shows  our 
discretized  model. 


M 


Figure  2.1:  Discretized  model  of  the  cantilever  beam 


where, 

=  Stiffness  of  the  ith  spring, 

—  Moment  corresponding  to  the  ith  actuator  (piezoceramic), 

U  =  Length  of  the  ith  link, 

(xi,  yi)  =  Position  of  the  center  of  mass  of  the  ith  link, 
o.i  —  Slope  of  the  ith  link. 

We  use  the  Euler-Lagrange  method  to  derive  the  equations  of  motion, 


d.8T. 

=  7t(ddJ 


dT  dU 

+ 


(2.1.1) 


dak  '  dak 

with  Mn+ 1  =  0.  The  potential  energy,  U,  and  kinetic  energy,  T,  are  obtained 
from  the  following  relations, 


(2.1.2) 


»=o  z 


(2.1.3) 
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Op  =  o,  QJJV+l  =  QljV, 


where, 

mi  =  Mass  of  the  ith  link, 

/,•  =  Moment  of  inertia  of  the  ith  link  with  respect  to  its  center  of  mass. 
Now,  we  compute  the  various  terms  in  the  Euler- Lagrange  equation, 


dT  ^  dxi  .  dxi  ^  ^  dyi  .  ,  %  ,  ,  . 

dcck  U  Uda>  dak  tk  U  dai  dak 


d  fdT  _  "  '  '  d2Xi  .  .  'dx;  dxj} 

dt  dak  ^  da^dai0^1  ^  da3  dak 

i  v^r  d2yi  .  .  v".  dyi  .  dyj , 

’  ,+§^' *'W 

4ws£“4&ci'1 

4m4£a4eStr]+IA’ 


(2.1.4) 


dU  dT 


dak  dak 


—  kk((Xk  &k— l)  kk+1  (o^+l  ®-k) 

4m4£d4^)] 


(2.1.5) 


Substituting  equations  (2.1.4)  and  (2.1.5)  into  the  Euler-Lagrange  equa¬ 
tion  (2.1.1  we  get, 


N 


M„-Mw  =  EK(  EE 


*'=i 


d2Xi  .  .  ^  dxi  ..  dxj 

i  dctjdcti^01  ^  da**3  dak 
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N 


d'yi 


-a; 


,=l  Ux  S'  daJda‘ 
-\-Ik&k  +  kk(ak  —  cvfc_i) 


(Xk)- 


u  dal 

kk+ 1  (^fc+1 


(2.1.6) 


the  position  of  the  center  of  mass  of  the  ith  link  can  be  obtained  from  the 
following  relations  (assuming  a  uniform  distribution  of  mass), 


i-i  j 

Xi  —  l{  COS  OLj  +  —li  COS  Cti 
i=i 

i-i  2 

Vi  =  ^2  *«'  s'm  ai  +  7)  k  sin  ai  ■ 
3= 1  Z 


(2-1.7) 

(2.1.8) 


To  simplify  the  equations  of  motion  (2.1.6)  we  need  the  partial  derivatives 
of  Xi  and  y,-  with  respect  to  ak  i.e., 


dxj 

dak 


0  i  <  k 

<  —  |  sin  ah  i  =  k 
—lk  sin  ak'  k  <  i 


dyi 

dak 


0  i  <  k 

i  ^  cos  ak  i  =  k 

—lk  cos  Ofc  k  <  i 


d2Xi 

da 2 


0 

—  ^  cos  aj 
—lj  cos  aj 


i  <  j 
i  =  j 
j  <  i 


(2.1.9) 


(2.1.10) 


(2.1.11) 
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0 


l  <  J 


d2yi 

da] 


—  ^  sin  aj  i  —  j  • 
— lj  sin  aj  j  <i 


(2.1.12) 


Now,  (2.1.6)  reduces  to, 


Mt-Mw  =  +  |1 

+hdk  +  h(oik  -  ak-i)  -  kk+i{ak+i  -  ak) 


N 


‘  1  d2X{  -j  dxi  ..  x  dxi 
'  daj^'dak1 


Er\ - N/'-/  Ust  ry  ..  v  ^ 

m*lE(  'a'l~2  ai  a"».  ai)  d'17' 


/-'v  da? 
j=i  uad 


Jfc— i 


v — \/d  ■„  dxk  ..  ~.dxk 

+miY,{—aj  +  —ai)^k 

^  d2x;  -2  ,  dxi  ..  dx; 

+  E  m4(-o~2  J 

i=fc+l  Wa‘ 


da;  1  dak 


+mk( 


d2xk  ;2  .  dxfc  ..  ,  dx* 


d«fc 


a£  + 


dak*k^dak 


rV^t/;  '2  ,  %  -  \  %l 

+  E  m4E(^-2 +  — ai)^J 


■  tj.i  da? 

t=K+l  J=1  3 


da ,  J'dafcJ 


k^/d2yk  -2  ,  dyk  dyk 

+m\§(d^^  +  dE^dE 

,  \/d2yi  *2  ,  %  ••  %1 

da?  da;  da*, 
,d2yfc  -2  ,  dyfc  ..  .  dy* 
+ra‘(M  ‘+fe 

+4<*fc  +  A;fc(afc  -  afc-i)  -  kk+i(ak+1  -  ak), 


(2.1.13) 


(2.1.14) 
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Mk  -  Mk~ i  = 


N  i- 1 

y  cosatjaf  -  lj  sin  ajdj)(-lk  sin  ak)] 

t=k+l  j= 1 

*-i  •  h 

—lk  YX- lj  cos  —  Ij  sin  ctjdj)(——  sin  ak) 

i=i  2 

N 


+  y  ^  cos  «.«i  ~  ^  sin a,at)(-4  sin  a*)] 


h 


1 


lk  ■•  \  /  lk  • 


+mk(-~-  cos  ^  sin  akdk)(—^  sin  afc) 

AT  »-l 

y  m,-[y  (— Zj  sin  aja]  +  lj  cos  aydy)(d-Zjk  cos  a*,)] 

»=A:+1  j=l 

*-i  .  •  h 

—lk  y  (— ■ lj  sin  oya|  +  lj  cos  ctjQ'j)(—  cos  a*) 
j=i  2 

w  ^  ■  b 

+  mi[(_77  s’n  aial  +  ^  cos  a> i&i)(h  COS  ajff)] 

*=fc+l  z  z 

+mk{-~  sin  dual  +  ^  cos  akak)(~  cos  a*.) 
d-/fc<*fc  +  Asjfc(o:jb  -  ttfc-i)  -  kk+i(ak+i  -  a*).  (2.1.15) 


After  some  cancellations,  the  above  equation  reduces  to, 


N  i- 1 

Mk  -  Mk+ 1  =  y  mi[y  ljlk{sin(ak  -  ocj)a )  +  cos(oijfc  -  a, )«.})] 
i=k+ 1  j— 1 

mk  k~1 

+— [y  ljh{sin(ak  -  aj)a]  +  cos(afc  -  ay)dy)] 

1  j=x 
JV  m. 

+  y  — l[lilk(sin(ak  -  a,)a-  -f  cos(afc  -  c*i)a,-)] 

i=rfc+l  Z 

r 2 

-(r  Ikcx.k  d-  kk{ctk  ctk—i )  —  ■  ^'k ) , 

jfc  =  1,  •••,  iV.  (2.1.16) 


These,  are  the  equations  of  motion  for  the  discretized  model. 
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In  order  to  apply  classical  control  strategies  we  will  have  to  linearize  them  about 
the  stationary  equilibrium  point  a.k  —  dk  —  0,  Jc  —  1 ,  •  •  ■ ,  N  . 

The  linearized  equations  are: 

Mt  -  Mt+ ,  =  e£m.  P'i  +  leh  Efc1  lM,  + \  Efat+i 

+  -  at-j)  -  fet+ifafc+i  -  a*), 

where  p  is  the  mass  density  per  unit  of  length  of  the  beam  and  we  made  the 
following  substitutions  for  the  moment  of  inertia,  h,  and  mass,  m^,  of  the  kth 
link  respectively, 


4  = 

mi  =  pli. 


2.2  An  Estimate  For  the  Stiffness  Constant  of 
the  kth  Link 

Up  to  this  point  we  derived  the  equations  of  motion  for  our  model.  The  next 
step  is  to  compute  the  various  constants  that  appear  in  (2.1.17)  .  One  of  these 
constants  is  the  stiffness  of  the  kth  link.  We  have  to  decide  on  some  criterion 
for  computing  these  constants.  One  reasonable  criterion  might  be  to  compute 
k{  such  that  the  deflection  of  the  tip  of  the  beam  at  the  steady  state  due  to  a 
force  F  applied  at  its  free  end  is  the  same  for  both  the  model  and  the  actual 
flexible  beam. 
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Figure  2.2:  Deflection  of  the  beam  due  to  a  force  F  applied  at  its  free  end 


The  deflection  of  a  cantilever  beam  due  to  a  force  F  at  its  free  end  is  [5], 


8  = 


FL3 

JeI' 


(2.2.1) 


where,  L  is  the  length  of  the  beam  and  E  is  the  modulus  of  elasticity.  For  our 
N-link  model  and  for  small  /5t-  in  the  static  situation  we  can  write, 


6  =  (2-2.2) 

i= 1  j= 1 

N 

kb  =  Fj^lr  (2.2.3) 

i=* 

Eliminating  from  the  above  two  equations  and  using  (2.2.1),  we  get, 


N  i  Nil 
i=l  j= 1  m=j  3 


V 3 


(2.2.4) 


We  also  assume  that  the  compliance  of  the  torsional  spring  corresponding  to 
each  link  is  proportional  to  its  length,  i.e., 


hli  —  c  i  —  1,  •  •  •  ,  N, 


(2.2.5) 
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where,  c  is  a  constant.  Combining  (2.2.4)  and  (2.2.5)  and  then  solving  for 
c  we  get, 


N  i  N  j  i  I . 

Lm 


EEE 


t=l  i=l  m=j 

3 El  N  i  N 


L3 

3  EV 


c  = 


L3  f 


EEEW4 


(2.2.6) 


t=l  J=1  771=7 

Alternatively,  one  may  consider  a  different  criterion  for  estimating  the  stiff¬ 
ness  coefficients:  find  ki  such  that  the  potential  energy  stored  in  the  model  is 
the  same  as  that  of  the  flexible  beam.  The  potential  energy  stored  in  torsional 
springs  is, 


u  =  Eh,  A2.  (2.2.7) 

i= 1  1 

The  potential  energy  stored  in  the  flexible  beam  is  obtained  from  the  follow¬ 
ing  equation: 

u=  2  L  Bl^)  dX-  (2'2'8) 

Again,  it  can  be  shown  that  the  transversal  deflection  of  the  flexible  beam 
due  to  a  force  F  applied  at  its  free  end  is  [5], 

u>(a;)  =  {2L3 -3L2x  +  x3)-^j  (2.2.9) 

where  x  is  the  distance  from  the  free  end.  Substituting  for  w( x)  from 
(2.2.9)  into  (2.2.8),  results: 

U  =  l/0iE/^[£(2L3-3L2x  +  I3)]2, 


F^L3 
6 El  ■ 
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Now  we  can  equate  (2.2.10)  and  (2.2.7)  to  get, 


F*L3  _  ST-h-S2 

6EI 


(2.2.10) 


Substitute  for  $  from  (2.2.3)  into  the  right  hand  side  of  the  above  equation: 


P  '2  J  3  i  N  P  N 

LJe  =  ±5>[f-V/i]2, 

*EI 


(2.2.11) 


and  finally  we  can  eliminate  fcj  from  the  above  equation  and  (2.2.5  to  obtain 
the  constant  c  , 


Thus, 


t  3  N  i  N 

w,  =  giW 

t— 1  J—X 


opjN  N 

c  =  (2.2.12) 

i=l  J=t 

Remark:  It  can  be  easily  shown  that  the  constant  c  goes  to  the  value  El 
as  the  number  of  rigid  links,  N  goes  to  infinity  while  the  length  L  remains 
constant. 


2.3  Relation  Between  Mi  and  the  Voltage 
Applied  to  the  Piezoceramic 

We  try  to  find  an  expression  for  Mi  in  terms  of  Vi  ,  the  voltage  that  applied 
to  the  ith  piezoceramic.  We  compute  Mi  such  that  at  the  static  situation  it 
results  in  the  same  amount  of  deflection  at  the  tip  of  each  link  as  the  amount  of 
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deflection  caused  by  applying  the  voltage  V  to  the  ith  segment  of  the  piezo 
material.  In  the  steady  state,  the  following  relation  holds  for  the  actual  flexible 
beam, 


g^EI  =  c,V  (2.3.1) 

where,  w  is  the  deflection  and  cp  is  a  constant  relating  the  voltage  applied  to 
the  piezo  material  to  the  moment  that  it  produces  [20].  Integrating  the  above 
equation  with  respect  to  x  two  times  we  get, 


(2.3.2) 


(2.3.3) 


Figure  2.3:  Deflection  resulting  from  the  moment  M{ 
where,  /ct-  has  already  been  computed. 
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2.4  Initial  State  Computation 


Suppose  our  model  for  the  cantilever  beam  is  subjected  to  a  force  F  at 
its  free  end.  In  order  to  perform  some  simulations  we  need  the  values  of 
cq,  i  =  1,  •  •  •  ,  TV  as  the  initial  state  for  the  dynamical  equations  of  motion. 
From  (2.2.1)  and  (2.2.3)  we  have: 


j?  N 

Pi  =  jT  XX'  ’ 

1  3- 1 


8  = 


FL3 

Jei ’ 


which  give  us  /?;  as  a  function  of  8  (deflection  at  the  tip), 


(2.4.1) 

(2.4.2) 


Pi  = 


1  3  EI8 
ki  L3 


N 


XX  i=  1,  •••  ,  N, 


and  ai  can  be  computed  using  the  following  relation, 


««■  =  i  Pi , 

_  3 EIS  1  y^/V  -j 

L3  1  2^j—m 


=  1, 


iV. 


(2.4.3) 


2.5  Simulation 

Based  on  equation  (2.1.17)  the  program  ’EqClMotion-m’  has  been  written  to 
derive  the  matrices  A  and  B  in  the  state  space  representation, 

x  =  Ax  +  Bu,  (2.5.1) 

which  we  need  for  simulation  and  control.  EqOfMotion.m  is  a  Mathematica 
macro.  As  input,  this  macro  accepts  n  =  number  of  links  and  physical  constants 
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of  the  beam  and  in  return  it  gives  the  numerical  matrices  A,  B ,  and  xq.  The 
output  of  the  program  is  stored  on  a  text  file  named  AandB.m  which  is  in  a 
format  that  is  readable  by  MATLAB.  Given  those  two  matrices  a  simulation  has 
been  performed  that  shows  the  free  oscillation  of  the  beam  when  it  is  released 
from  the  initial  rest  with  8  —  4  cm. 

The  results  are  showm  on  figure  2.4  and  figure  2.5  for  N  =  3  and  N  = 
4  respectively.  The  upper  graphs  are  the  slopes  of  each  rigid  link  versus  time 
and  the  lower  ones  shows  the  time  derivative  of  the  slope  of  the  rigid  links.  Each 
curve  corresponds  to  one  link. 


Figure  2.4:  Free  oscillation  of  the  beam  N  =  3 

One  may  compare  the  lowest  mode  of  oscillation  of  the  Galerkin  model  ob¬ 
tained  above  (which  is  easily  readable  from  the  deflection  curves)  and  compare 
it  with  the  one  that  can  be  obtained  analytically  by  solving  the  Euler- Bernoulli 
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Figure  2.5:  Free  oscillation  of  the  beam  N  =  4 

equation  to  verify  the  correctness  of  the  model  and  simulation.  In  the  next 
chapter  we  follow  a  different  algorithm  to  get  a  finite  dimensional  model  for  the 
flexible  beam. 
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Chapter  3 


Matching  Pursuit  with  Rational 
Wavelets  for  Nonparametric 
Estimation  of  Stable  Systems 


The  matching  pursuit  algorithm  developed  by  S.  Mall  at  and  Z.  Zhang  [18] 
has  been  exploited  to  obtain  a  rational  transfer  function  approximation  to  a 

set  of  empirical  data  obtained  from  say,  an  experiment,  or  even  to  reduce  the 

> 

order  of  a  transfer  function  which  is  real-rational  already.  In  this  chapter  we 
follow  closely  the  work  done  by  Y.  C.  Pati  in  his  Ph.D.  thesis  [15],  which 
basically  addresses  the  same  problem  but  with  a  different  solution.  Using  this 
method,  a  priori  knowledge  of  time-frequency  localization  properties  of  the  given 
transfer  function  enables  us  to  obtain  a  “good”  model  for  the  physical  system 
under  consideration.  One  advantage  of  this  algorithm  is  its  simplicity  which 
in  turn  makes  it  fast  (e.g.  compared  to  the  least  squares  method  discussed  in 
[15]).  As  a  practical  application  of  this,  we  considered  the  problem  of  finding  a 
finite  dimensional  approximation  to  the  transfer  function  of  flexible  beam  with 
°This  work  has  been  done  under  the  supervision  of  Dr.  Y.  C.  Pati. 
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piezoceramic  actuators.  The  performance  of  the  matching  pursuit  algorithm  is 
compared  to  the  Laguerre  method  and  also  the  one  developed  in  [15]. 

3.1  Introduction 

The  nonparametric  estimation  problem  has  been  addressed  in  many  ways  and 
different  methods  have  been  devised  to  solve  the  problem.  Obviously  one  always 
seeks  a  low  order  approximation  while  achieving  a  certain  degree  of  accuracy. 
It  is  also  desirable  to  obtain  an  approximation  that  can  be  easily  realized.  We 
use  the  algorithm  developed  by  S.  Mall  at  and  Z.  Zhang  [18]  that  decomposes 
any  signal  into  a  linear  combination  of  waveforms  that  belong  to  a  redundant 
dictionary  of  functions. 

Clearly,  the  choice  of  the  dictionary  has  a  great  deal  to  do  with  how  compact 
the  decomposition  is  going  to  be.  For  the  system  identification  problem  one  of 
the  requirements  is  the  real-rationality  of  the  approximation  (a  real-rational 
function  is  a  rational  one  with  real  coefficients).  At  least  this  is  important  when 
we  want  to  use  the  method  to  represent  an  infinite  dimensional  system  with  a 
finite  dimensional  one. 

We  also  may  use  our  a  priori  knowledge  of  time-frequency  localization  prop¬ 
erties  of  the  signal  to  select  the  elements  of  the  dictionary  in  an  efficient  fashion. 

Motivated  by  these  two,  rational  wavelets  seem  to  be  a  good  choice  for  the 
purpose  of  decomposition  of  a  transfer  function  because: 

•  1.  They  are  rational. 

•  2.  They  capture  the  localization  properties  of  the  signal. 
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•  3.  It  is  easy  to  realize  a  wavelet  decomposition  of  the  signal  (elements  of 
the  dictionary  are  just  translation  and  dilation  of  the  mother  wavelet.)). 

Of  course,  the  fact  that  the  approximation  given  by  this  method  is  real-rational 
is  something  that  must  be  proven  and  we  will  show  that  in  section  3. 

One  thing  that  should  be  pointed  out  is  that  although  a  poor  choice  will 
affect  the  compactness  and  the  degree  of  accuracy,  the  elements  of  the  dictio¬ 
nary  can  be  virtually  any  set  of  functions  and  there  are  no  constraints  such  as 
orthogonality  on  them.  The  only  requirement  is  that  they  must  be  unit  normed. 

3.2  Preliminaries 

Before  introducing  the  algorithm,  we  need  to  develop  some  notations  which  are 
almost  standard  but  yet,  for  the  sake  of  completeness,  worth  mentioning.  The 
inner  product  of  two  signals  /,  and  g  in  L2(f?)  is  defined  to  be 


/+oo 

f{t)g(t)dt, 

-OO 

which  leads  to  the  usual  L2  norm, 

/+CO 

\\m\\2dt. 

-oo 

If  $(cu)  is  the  analyzing  wavelet  $m,n(s)  is  defined  by, 

$m,7i(s)  =  a™/2$(a™s  -  inbo),  (3.2.1) 

where  a0  and  b0  are  dilation  and  translation  step  sizes,  respectively  and  m  and 
n  are  called  dilation  and  translation  levels.  Note  that  it  is  necessary  to  put  the 
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imaginary  number  i  in  (3.2.1)  for  n  to  be  the  translation  index.  To  see  this,  let 
us  substitute  iu>  for  s  in  (3.2.1), 

^m,n(iw)  =  a™l2§{a™iio  -  inb0) 

=  -  aQmnb0)). 

3.3  The  Algorithm 

First,  we  describe  the  algorithm  in  its  general  form  and  its  application  to  rational 
transfer  functions  follows  as  a  special  case.  In  the  initialization  part  of  the 
algorithm,  a  set  of  translation  and  dilation  levels  has  to  be  chosen  that  is  referred 
to  as  dictionary.  This  is  the  set  from  which  the  algorithm  select  the  translation 
and  dilation  levels.  The  apriori  knowledge  of  time-frequency  localization  of  the 
signal  should  help  to  choose  this  set.  Here  the  objective  is  to  approximate  /, 
an  element  of  the  Hilbert  space,  by  a  linear  combination  of  elements  of  the 
dictionary. 

•  Step  0:  Select  the  analyzing  wavelet  $  with  norm  equal  to  one. 

•  Step  1:  Choose  the  finite  set  D  of  translation  and  dilation  levels. 

•  Step  2:  Set 

R°f  =  /;  i  =  I- 

•  Step  3:  Find  indices  (m;,ti;)  £  D  such  that, 

l(R ‘-'A  <W,)I  >  l( (3.3.1) 

V(m,  n)  €  D. 
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•  Step  4-'  Compute  the  approximation  at  level  i  via, 


(3.3.2) 


3.4  The  Theory  Behind  The  Algorithm 

First  note  that  the  residuals  satisfy  the  following  recursive  equation, 

fl"/  =  +  R"+'f.  (3.4.1) 

where,  $7n  is  the  selected  wavelet  at  the  nth  iteration.  Evaluating  the  inner 
product  on  both  sides  with  $7n  and  noting  that  $7n  is  normed  one,  proves  that 
Rn+1  f  is  orthogonal  to  $7n  which  in  turn  enables  us  to  write, 

IIK/ll2  -  I (Rnf,  <M|2  +  ||i?”+1/||2.  (3.4.2) 

This  shows  that  the  residual  Rnf  is  monotonically  decreasing.  In  fact  the 
following  theorem  states  that  the  series: 

OO 

/  =  £(/?"/,  47„)4.7„  (3.4.3) 

n=0 

converges  to  the  projection  of  /  into  the  subspace  spanned  by  the  dictionary. 


/«'  —  fi  ^mk,nk)$mk,nk 

k= 0 

If  satisfactory,  stop.  Otherwise  continue. 

•  Step  5:  Compute, 

R'f  =  R^f-iR^f, 

increment  i,  and  go  to  step  3. 
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Theorem  [18]:  Let  H  be  a  Hilbert  space  and  /  G  H.  The  residual  Rm  f  defined 
by  the  recursion  (3.4.1)  satisfies 


lim  \\Rm f  —  Pw/||  —  0.  (3.4.4) 

m— +oo 

Hence 


p vf  =  £{*“/,  (3-4.5) 

71—0 

and 


IIPv/ll2  =  £1(^7,  (3.4.6) 

n= 0 

Where,  V  is  the  space  spanned  by  the  vectors  in  D,  W  is  its  orthogonal  com¬ 
plement  in  H,  and  Py  and  P w  are  orthogonal  projections  of  /  into  V  and  W 
respectively.  When  V  =H,  we  have  Py/  =  /  and  the  matching  pursuit  recovers 
/• 

But,  if  we  want  to  use  this  method  to  approximate  a  stable  transfer  function 
F  with  real  coefficients  there  is  one  more  thing  that  needs  to  be  shown.  That 
is,  we  have  to  show  that  the  projection  of  F  into  the  space  V  is  a  real  rational 
function.  Formally,  the  space  of  transfer  functions  of  stable  systems  is  the  Hardy 
space  H2(n+))  where  n+  is  the  right  half  plane  .  Also,  RH2(n+)  is  a  subset  of 
H2(n+))  that  consists  of  those  elements  that  are  rational  functions  in  s  with 
real  coefficients.  So,  we  need  to  prove  that  Py/  =  /  GRH2(n+)  . 

Rationality  is  not  a  problem  because  the  analyzing  wavelet  is  rational  to 
begin  with  (we  choose  it  to  be  this  way).  But  note  that  $m>n(5)  =  a£/2§(a™s  — 
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inb0)  is  not  the  Laplace  transform  of  a  real  valued  function  (unless  n  —  0).  It 
can  be  easily  shown  that  [15]  Gm,n(s )  defined  as, 

Gm’n(s)  =  a$m,»(s)  +  a$m,-«(s) 

is  a  real  rational  function  in  H2(IT+). 

Thus,  to  guarantee  the  real  rationality,  for  every  wavelet  in  the  space 

V ,  corresponding  wavelet  with  the  translation  index  —n  should  be  present  in  V 
as  well.  In  fact  by  doing  this,  the  matching  pursuit  process  picks  the  coefficients 
corresponding  to  wavelets  4>min  and  complex  conjugate  of  each  other  to 

make  P vf  real  rational. 

Hereafter,  we  use  $(.)  instead  of  <h(L)  to  simplify  the  notation.  Recall  that 
if  the  inverse  Fourier  transform  of  F(u>)  is  a  real  valued  function  then, 

F(—u))  =  FH-  (3.4.7) 

We  will  use  this  in  the  proof  of  the  following  proposition. 

Proposition  1  The  following  relation  holds, 

(F,  $TO)_n)  =  (F,  (3.4.8) 

Proof: 


(f,  =  I^F^)a^i(a^w-nb„) 

=  F(»)aZ'%(nk,  - 

=  {F  *  $m,o){nb0/aZ) 

—  •^r[/(^)^m,o(^)]L=n6oa~m- 

Similarly, 
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[f, 


(3,4,9) 


Since  f(t )  and  are  both  real  we  have, 

and  this  proves  the  statement  of  the  proposition.  □ 

The  next  proposition  is  what  we  need  to  justify  that  the  approximation  given 
by  the  matching  pursuit  algorithm  converges  to  a  real  rational  function. 

Proposition  2  Let  V  be  such  that  3>TO)„  £  V  implies  that  $m-n  £  V  for  all 
m,  n  £  ,D,  then  Py /  is  real  rational. 

Proof:  Since  Py /  is  the  orthogonal  projection  of  /  into  V  we  have, 

(/  -  P v/,  5)  =  0,  V#  £  V. 

In  particular  the  following  equalities  hold: 


(/(w),  $m,n(w)  +  =  (Pv/H,  $m,n(w)  +  $m,-n(w))  (3.4.10) 

(/(w),  $m,»M  ~  $TO|-n(t*>))  =  (Py/(w),  $TO,„(a;)  -  (3.4.11) 

From  (3.4.7)  we  see  that  <hTOjn(a))  +  $m-n(u)  is  a  real  rational  function  and 
by  (3.4.7)  we  can  write, 

®m,7i(— u^)  +  <hTO)_„(— w)  =  5>m)n(u>)  -f  Omj_n(o;).  (3.4.12) 
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Replace  uo  by  —u  in  (3.4.11)  and  use  (3.4.12)  and  the  fact  that  f (co)  has  a 
real  valued  weighting  pattern  (therefore  /(— u>)  —  /(co)),  to  get 

(Py/(-w),  $m,„(-w)  +  $m,-»(-w))  =  (/(-w),  $m,n(-w)  +  $m,-n(-w)) 

=  (/(w)?  ^m,n(w)  4~  ^m,-n(w)) 

(3.4.13) 

—  5  4“  n(^)) 

=  (Pv/(w)>  ^*TO,Ti(w)  +  ^m,-n(^))- 

Therefore, 

(Py/(-W)-W7R,  *m,nM  +  =  0,  (3.4.14) 

or,  after  taking  complex  conjugate, 

(Py/(-u)-Py/(u),  $m,n(^)'f  $m,-n(w))  =  0  (3.4.15) 

V  (m,n)  G  D. 

It  can  be  shown  that  the  following  holds: 

k*)  n{  ^)  =  ^m,n(w)  4"  n(w)*  (3.4.16) 

Starting  with  (3.4.11)  and  taking  the  same  steps  as  we  took  for  (3.4.11),  we 
obtain 

(Py/(— oO-Py/M,  $m,B(u)  -  =  0,  (3.4.17) 

V  (m,n)  6  D. 
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Add  equations  (3.4.16)  and  (3.4.18)  up  to  get 


(P vf(-co)  -  Py/M,  =  0,  (3.4.18) 

V  (m,n)  G  D. 

Equations  (3.4.16)  and  (3.4.18)  can  also  be  rewritten  as 

(P v/(-w)  -  Pv/H,  $m,n(-w)  +  $m,_n(-w))  =  o,  (3.4.19) 

V  (m,  n)  G  D, 

(MH  -  Pv/H,  $TOl«(-w)-$TO,-n(-w)>  =  0,  (3.4.20) 

V  (m,n)  G  D, 

which  are  obtained  by  taking  the  complex  conjugate  once  and  using  real  ra¬ 
tionality  of  $m,„(w)  +  $TO)_n(a;)  and  equation  (3.4.16).  Again,  add  equations 
(3.4.20)  and  (3.4.21)  up  to  obtain 


(P vf(-to)  -  Pv/M,  $m,n(-u;))  =  0,  (3.4.21) 

V  (m,  n)  G  D. 

Notice  that  P yf  and  Py/(-w)  can  be  expressed  as  a  linear  combination  of 
vectors  in  V.  i.e. 


pv/H  = 

Pv/(-w)  =  Eue/$y(^)  I C  D. 

From  equations  (3.4.22),  (3.4.19),  and  (3.4.22)  we  conclude  that, 


(3.4.22) 
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(P vf(-uj)  -  P y/(aj),  P vf(-U>)  -  P v/(w))  =  o. 


(3.4.23) 


Hence, 


Pv/(-w)  —  Py/(w)  almost  everywhere.  (3.4.24) 

which  together  with  continuity  of  Py/(u>)  results  that  the  equality  holds  ev¬ 
erywhere.  This  is  the  necessary  and  sufficient  condition  for  Py/  having  a  real¬ 
valued  weighting  pattern.  □ 

3.5  Computational  Aspects 

There  are  a  few  points  worth  noting  that  make  the  algorithm  more  efficient  as 
far  as  storage  and  speed  is  concerned.  First  note  that  in  order  to  choose  the 
indices  (m,-,  n,)  in  the  third  step,  and  also  to  compute  the  approximation  at  step 
4,  we  only  need  the  inner  product  of  the  residuals  with  the  wavelets  (and  not 
the  residuals).  Thus,  in  step  5  one  may  directly  compute  the  inner  products 
instead,  i.e. 

(i?7,  *m,»>  =  (3.5.1) 

Note  that  by  doing  this  we  actually  do  not  have  to  do  any  integration  for 
computing  the  inner  products  (other  than  just  one  set  of  integration  to  get 
(/,  $m,n)).  As  far  as  storage  is  concerned,  note  that  at  each  step  we  only  have 
to  keep  one  set  of  inner  products  of  the  residuals  with  the  wavelets  and  as  we 
proceed,  we  can  throw  away  the  previous  ones. 
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The  third  point  is  that  the  correlation  matrix  (with  the  elements  ($k,h 
$mtn))  can  be  computed  a  priori  and  in  fact  this  can  be  done  analytically: 


$ma,n2M)  =  /  $mi ,ni (u>)a™2/2$(cua™2  -  n2b0)duj 
=  <h2/2  f  §muni(u))$(n2bo  -  uj™2)(1lo 


—  ( ,n\  *  ^m2,n2)|a;_Ii2^a 

w —  mo 

“o 

=  •^~{4>Tni,ni  (f)^Tn2,o(f  ))  |w=H2^a  • 

Qo  2 

Also,  4>m,n(t)  can  be  written  in  terms  of  <f>(t)  as, 


,  /.x  _  -m/2  ,(  t  x 

<^m,n(0  —  a0  e  0  m)‘ 

a0 


Substitute  for  from  (3.5.3)  into  (3.5.2)  to  get, 


(^mi.ni  (k-*),  *^7712,712  (^-O)  —  a\ 


t 


(3.5.2) 


(3.5.3) 


F(e-a°  m2))\u=zga._z£o. (3.5.4) 

uo  uo  »n  2  °n  1 


3.6  Results 

The  matching  pursuit  scheme  has  been  used  to  obtain  the  approximation  to 
several  frequency  responses  and  in  each  case  we  compared  its  performance  with 
the  Laguerre  method  and  wavelet  decomposition  using  least  square  [15].  As  an 
example  consider  the  second  order  system  with  delay  which  has  been  examined 
in  [15].  The  transfer  function  for  this  system  is 

4.94e~2s 

s2  +  1.25s  +  0.406 

Figure  (3.1)  shows  the  magnitude  of  the  wavelet  coefficients  obtained  from 
the  matching  pursuit  algorithm.  The  analyzing  wavelet  is  again  given  in  Eq.(6.8.1) 
As  it  can  be  seen  in  the  figure,  very  few  coefficients  have  a  significant  magnitude 
and  this  is  a  measure  of  how  well  the  frequency  response/impulse  response  is 
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localized  in  time-frequency  domain.  This  localization  property,  in  turn,  enables 
us  to  get  a  fairly  low  degree  for  the  approximation. 


15  80 


Translations 


Dilations 


Figure  3.1:  Magnitude  of  the  wavelet  coefficients  (delay  system) 

The  approximation  results  together  with  the  true  values  of  the  frequency 
response  are  shown  in  figure  (3.2).  The  model  order  is  16  here  and  the  number 
of  wavelets  that  has  been  used  is  only  5  (each  wavelet  with  nonzero  translation 
increase  the  order  by  four  and  those  with  zero  translation  only  increase  the 
model  order  by  two). 

For  the  purpose  of  comparison,  the  normalized  error  (in  time  domain)  versus 
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Real  part  Imaginary  part 


Figure  3.2:  Approximation  to  the  delay  system  with  model  order  1G 

the  model  order  is  plotted  in  figure  (3.3).  It  can  be  seen  from  figure  (3.3)  that 
the  matching  pursuit  algorithm  outperforms  both  methods  in  a  relatively  wide 
range  of  model  order  and  it  is  only  with  high  order  models  that  the  least  square 
method  beats  the  matching  pursuit.  In  this  figure,  the  Laguerre  data  is  taken 
from  [15] 

In  figure  (3.4)  is  again  the  magnitude  of  the  wavelet  coefficients  for  the 
cochlear  filter  (it  is  a  filter  in  the  human  ear  that  separates  different  frequencies 
in  the  receiving  wave).  The  same  set  of  graphs  are  shown  for  in  figure  (3.5).  As 
before,  the  time-frequency  localization  properties  are  clear  from  figure  (3.4). 
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Figure  3.3:  Normalized  time  domain  approximation  error  versus  model  order  for 
different  algorithms  (delay  system)  (Laguerre  data  from  [14]) 

3.7  Application 

As  a  practical  application  we  consider  the  problem  of  approximation  of  the 
transfer  function  relating  the  input  voltage  to  the  actuating  piezo  ceramic  to 
the  output  voltage  of  the  sensors  (which  are  piezo  ceramics  again).  The  setup 
for  this  experiment  is  shown  in  figure  (6.7  ). 

The  empirical  transfer  function  has  been  obtained  through  simulation  and 
application  of  spectral  analysis. 
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Figure  3.4:  Magnitudes  of  wavelet  coefficients  (cochlear  filter) 

Due  to  the  presence  of  resonances  the  frequency  response  turns  out  to  be 
fairly  well  localized  in  the  frequency  domain,  which  is  a  desirable  factor  when 
we  want  to  apply  the  matching  pursuit  algorithm. 

Figure  (3.7)  shows  the  result  of  an  approximation  of  degree  28.  The  dashed 
line  represents  the  empirical  frequency  response,  and  the  solid  line  is  the  the  one 
obtained  by  the  degree  28  approximation.  Note  that  the  governing  model  for 
the  problem  is  actually  a  PDE  and  thus  it  is  not  a  finite  dimensional  system. 
The  analyzing  wavelet  is  again  the  one  given  in  Eq.(6.8.1). 

It  is  also  interesting  to  get  a  high  order  approximation  and  observe  that 
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Real  part  Imaginary  part 


Figure  3.5:  Approximation  to  the  cochlear  filter  with  model  order  16. 

only  a  few  of  the  wavelet  coefficients  are  of  significant  magnitude.  Figure  (3.8) 
reflects  this  fact. 

3.8  The  Experiment 

In  addition  to  previous  applications  we  arranged  an  experimental  setup  consist¬ 
ing  of  a  flexible  cantilever  Aluminum  beam,  piezo  material,  spectrum  analyzer, 
and  a  power  amplifier  together  with  two  power  supplies.  The  objective  here, 
as  in  its  simulated  version,  is  to  find  a  finite  dimensional  model  that  approx- 
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Figure  3.6:  Normalized  time  domain  approximation  error  versus  time  for  differ¬ 
ent  algorithms  (cochlear  filter)  (Laguerre+Least  square  from  [14]). 

imates  the  transfer  function  from  input  voltage  to  the  power  amplifier  to  the 
output  voltage  measured  by  the  piezo  ceramics.  A  schematic  diagram  of  the 
experimental  setup  is  shown  in  fig. (3.9)  . 

The  spectrum  analyzer  (HP  3566A)  [6]  applies  a  swept  sine  wave  to  the  input 
of  the  power  amplifier  and  senses  the  voltage  generated  by  the  sensor.  The  gain 
of  the  power  amplifier  is  about  43  and  it  has  a  fairly  large  bandwidth  [3].  The 
frequency  response  obtained  in  this  way  is  actually  the  one  from  the  input  of  the 
amplifier  to  the  sensor  voltage.  The  output  voltage  of  the  spectrum  analyzer  is 
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Real  part  Imaginary  part 


Figure  3.7:  Approximation  of  order  28  to  the  beam  transfer  function 

adjustable  and  has  been  set  to  get  an  output  voltage  of  about  120  volts  pk  —  pk 
at  the  output  of  the  amplifier.  At  frequencies  below  80 Hz  the  output  is  distorted 
and  it  is  so  small  that  the  60 Hz  noise  of  the  line  is  more  significant,  so  we  start 
from  80 Hz.  The  Bode  plot  obtained  from  the  experiment  is  shown  in  fig. (3. 11). 
The  peaks  in  the  magnitude  plot  represent  the  resonance  frequencies. 

Having  the  empirical  data  at  hand,  the  next  step  is  to  feed  it  to  the  program 
and  get  the  approximation.  Here  there  are  a  few  parameters  that  we  should 
choose  before  running  the  program.  First,  we  must  choose  the  translation  step- 
size  6q.  There  is  a  limit  on  how  large  we  can  choose  this  parameter  and  it  is 


37 


Wavelet  coeff.  magnitude 


2 

1.5 

1 

0.5 

0 

-0.5 

1 


0 


Figure  3.8:  Magnitudes  of  wavelet  coefficients  (experimental  set  up) 

determined  by  a  limit  beyond  which  the  set  of  wavelets  cease  to  be  a  frame  (see 
[15]  for  a  definition  of  frames).  As  the  analyzing  wavelet  again  we  chose  the 
second  order  system, 

s)  = _ 1 _ 

1  o+7)2+r 

with  the  parameters  as  in  Table  (3.1).  The  upper  limit  for  b0  to  preserve  the 
frame  property  is  in  this  case  is  about  17  [15].  In  Table  (3.2)  we  have  listed 
the  selected  elements  of  the  dictionary  for  an  approximation  of  order  46.  Note 
that  each  wavelet  with  nonzero  translation  adds  up  four  to  the  order  (because 
we  have  to  include  the  negative  translation  too)  and  those  with  zero  translation 
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7 

o0 

bo 

1.0 

5.0 

2.0 

16.37 

Table  3.1:  Parameters  of  the  analyzing  wavelet  and  translation-dilation  step 


Figure  3.9:  Schematic  diagram  of  the  experiment  set  up 
just  add  two.  The  approximation  has  the  following  form: 

M  =  En#o  J2m(Cmnaf2g(af2s  -  inb0)  +  C^a^/2g(a^2s  +  inbo))+  g 
Em  Cmoao/2g{a™/2s). 

The  frequency  response  of  the  approximation  corresponding  to  Table(3.2), 
along  with  the  empirical  frequency  response  is  plotted  in  Fig. (3. 12).  Fig. (3. 14) 
shows  the  magnitude  of  the  wavelet  coefficients  for  a  high  order  approximation. 

3.9  Conclusion 

The  matching  pursuit  method  is  a  simple  yet  efficient  algorithm  to  compactly 
represent  a  waveform/signal.  Here,  with  a  little  modification,  we  applied  this 
method  to  the  problem  of  obtaining  an  approximation  to  a  given  frequency  re¬ 
sponse.  Through  comparison,  it  has  been  shown  that  this  method  out  performs 
the  classical  Laguerre  method  for  system  identification  and  gives  a  lower  nor¬ 
malized  error  for  the  same  degree  of  the  approximation.  There  is  a  natural 
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Dilation 

Translation 

Cm  7i 

level 

level 

i — » 

59 

-0.69032973051  -  3.84848356247^ 

0 

52 

3.12591791153  +  2.19522929192i 

9 

0.40908780694  +  3.03336858749* 

29 

2.59771537781  +  0.25139620900i 

-1 

15 

-3.00720357895  -  7.58607721329; 

-1 

7 

2.81444954872  -  0.27529984713i 

-1 

20 

1.73732197285  —2.21613240242* 

-2 

13 

4.71194553375  +  23.73843002319* 

-3 

5 

1.66028666496  -  9.23096656799* 

-4 

3 

-3.34253716469  -  0.5581 1583996i 

-4 

2 

0.47969451547  -  6.92282009125* 

-5 

2 

5.80866384506  +  0.55891269445* 

-6 

1 

-2.26374673843  -  2.98484420776* 

Table  3.2:  Selected  elements  of  the  dictionary  and  their  coefficients 
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Rf  =  1 0Okfi 


Figure  3.10:  Piezo-electric  transducer  drive  circuit 

question  here  that  one  may  ask  about  this  method,  i.e.  since  the  approximation 
actually  converges  to  the  projection  of  the  signal  into  the  space  spanned  by  the 
waveforms  in  the  dictionary,  why  don’t  we  simply  compute  the  projection  using 
pseudo  inverse  and  in  this  way  get  the  answer  in  one  shot?  There  are  two  rea¬ 
sons  that  we  do  not  want  to  do  this;  Firstly,- computation  of  the  pseudo  inverse 
is  very  expensive  and  requires  lots  of  operations  specially  when  we  are  dealing 
with  extremely  redundant  dictionary.  The  second  reason,  which  is  really  more 
important,  is  that  our  objective  is  to  pick  those  waveforms  in  the  dictionary 
that  are  more  important  in  the  sense  that  they  can  represent  the  structure  of 
the  signal  best.  In  this  way  we  have  taken  advantage  of  the  redundancy  of  the 
set  of  waveforms  and  at  the  same  time  if  we  want  to  only  choose,  say  n  number 
of  waveforms  to  represent  the  signal  we  know  when  to  stop  where  as  in  the  case 
of  pseudo  inverse  there  is  no  clear  way  to  do  this. 


42 


44 


Real  part 


Modulus 


250 

200 

150 

100 

50 

0 


. Approx. 

- True 

1 

i 

:  1 

| 

!  j  ,  a  A _..J 

i 

_ 

0 


500 


1000 


1500 


T?: _ o  io-  tt:  _ 


mors 


moo: 


Imaginary  part 


Dilations 


Translations 


Figure  3.14:  Magnitude  of  the  wavelet  coefficients  corresponding  to  the  experi¬ 
mental  transfer  function 


Chapter  4 


Design  of  Controller  for  the 
Constrained  System 


Like  most  of  the  systems  in  the  real  world,  there  is  a  limit  on  the  amplitude 
of  the  input  that  can  be  safely  applied  to  the  actuators.  In  the  case  of  flexible 
beam,  these  limits  are  imposed  partly  by  the  saturation  limit  on  power  amplifiers 
driving  the  piezoceramics  and  partly  by  the  break  down  voltage  of  the  piezo  ma¬ 
terial  itself.  One  approach  to  controller  design  would  be  to  define  a  cost  function 
and  solve  a  constrained  optimization  problem  to  compute  the  optimum  input. 
Although  this  gives  the  “optimal”  solution,  it  has  the  disadvantage  that  both 
computation  and  storage  of  the  optimum  input  is  expensive.  As  an  alternative 
procedure  we  consider  the  method  proposed  by  Gutman  and  Hagander  [13]  and 
formulate  the  vibration  suppression  problem  into  a  form  suitable  for  application 
of  this  control  strategy.  The  results  of  this  method  are  then  compared  with  the 
LQR  method. 
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4.1  Application  of  the  LQR  method 


As  mentioned  before,  we  assume  that  there  is  a  constraint  on  the  magnitude  of 
the  maximum  and  minimum  voltages  that  can  be  applied  to  the  input  of  the 
power  amplifiers  driving  the  piezoceramic  actuators.  As  a  basis  for  comparison 
first  we  use  the  LQR  method  for  obtaining  the  suitable  stabilizing  feedback  gain 
L.  The  penalty  on  inputs  has  been  increased  until  the  magnitude  of  inputs 
satisfy  the  constraints  at  all  times.  The  results  of  the  simulation  are  shown  on 
figure  4.1.  Here,  we  used  the  linearized  Galerkin  model  obtained  in  chapter  2. 
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Figure  4.1:  Controlled  vibration  of  the  beam  (Linear  Quadratic  Regulator) 

Shown  in  figure  4.2  and  figure  4.3  are  the  results  of  the  simulation  for  the 
cases  that  a  disturbance  torque  in  the  form  of  a  sinusoid  at  the  frequency  of  the 
largest  and  the  smallest  modes  is  applied  at  the  first  joint,  respectively. 
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The  effect  of  disturbance 


Figure  4.2:  Controlled  vibration  of  the  beam  under  the  disturbance  torque  (Lin¬ 
ear  Quadratic  Regulator) 

Now  let  us  put  a  saturation  block  at  the  input  terminal  of  the  system  and 
observe  the  effects  of  saturation  with  the  LQR  as  the  controller  as  we  increase 
the  amount  of  the  initial  deflection  of  the  free  end  of  the  beam.  Simulation 
results  are  shown  in  figure  4.4. 

4.2  The  Gutman-Hagander  Method 

In  this  section,  we  use  a  control  strategy  devised  by  Gutman  and  Hagander  [13] 
for  controlling  systems  subject  to  input  constraints.  Although  this  method  is 
not  optimal,  it  guarantees  the  asymptotic  stability.  Hereafter,  in  this  chapter 
we  will  use  the  same  notations  and  symbols  as  in  [13]  and  the  reader  is  referred 
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The  effect  of  disturbance 


Figure  4.3:  Controlled  vibration  of  the  beam  under  the  disturbance  torque  (Lin¬ 
ear  Quadratic  Regulator) 

to  that  paper  for  a  detailed  description  of  the  method  and  definition  of  various 
symbols  that  we  are  going  to  use.  The  general  idea  is  to  first  stabilize  the 
system  by  a  low-gain  linear  state  feedback.  Then,  quadratic  Lyapunov  function 
is  found,  on  the  basis  of  which  another  linear  state  feedback  is  computed.  The 
two  controls  are  added  and  passed  through  the  saturation  element  Figure  (4.5) 
shows  the  block  diagram  of  such  a  controller. 

The  algorithm  consists  of  five  steps: 

•  Step  1:  Determine  the  set  of  initial  conditions,  D.  This  is  the  set  from 
which  we  want  the  controller  steer  the  state  of  the  system  to  zero.  One 
obtains  an  estimate  for  this  set  by  a  knowledge  of  the  physics  of  the  system. 
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Input  voltage  [voltsj  Deflection  [Deg] 


LQR  with  saturation 


Figure  4.4:  The  effect  of  saturation  on  the  LQR  controller 

•  Step  2:  Find  the  low-gain  stabilizing  feedback  gain  by  solving  an  LQR 
problem.  Increase  the  control  penalty  until  the  inputs  satisfy  the  con¬ 
straint  for  x  starting  in  D. 

•  Step  3:  Find  the  positive  definite  matrix  P  by  solving  the  following  Lya¬ 
punov  equation, 


PAC  +  A?P  =  -Q,  (4.2.1) 

where,  Q  is  a  positive  definite  matrix  that  can  be  considered  as  a  control 
parameter  and  Ac  —  A  +  BLT.  Since  Ac  is  stable  this  equation  has  a 
unique  positive  definite  solution. 

•  Step  p.  Let  E  be  the  set  of  state  variables  such  that  the  control  u  =  LTx 
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Figure  4.5:  Block  diagram  of  a  Gutman-Hagander  type  controller. 


satisfies  the  constraint.  Check  to  see  if  the  following  inequality  holds, 

sup  xTPx  <  min  xT Px.  (4.2.2) 

If  not  either  go  back  to  step  3  and  choose  another  matrix  Q,  or  change  the 
feedback  gain  L,  or  choose  a  more  restrictive  set  of  initial  conditions  D. 

•  Step  5:  Set  up  the  control  u  and  tune  the  parameter  /F,  which  is  a  diagonal 
positive  definite  matrix, 

u  =  sat[(LT  -  KBtP)x).  (4.2.3) 

In  the  following  subsections,  we  formulate  the  problem  of  vibration  damping 
such  that  fits  the  Gutman-Hagander  (G-H)  algorithm. 

Identifying  the  Initial  Condition  Set  D 

Assume  that  the  objective  is  to  damp  the  vibration  of  the  beam  resulting  from 
deflecting  the  tip  of  the  beam  and  then  releasing  it.  Simulation  results  of  free 
oscillation  can  be  used  to  obtain  the  region  from  which  we  want  to  steer  the 
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state  to  zero.  i.e.  by  looking  at  the  maximum  and  minimum  values  that  each 
state  variable  takes  during  free  oscillation,  we  can  get  an  estimate  of  the  range 
of  values  on  which  state  variables  take  values.  This,  of  course,  gives  a  conserva¬ 
tive  and  rather  large  estimate  for  the  set  D  because  the  extremum  of  all  state 
variables  do  not  occur  at  the  same  instant  of  time.  Thus,  we  get  inequalities 
like  the  following, 


rrii  <  Xi  <  Mi  i=  1,  •••  ,  N.  (4.2.4) 

Constraints  on  Inputs 

As  mentioned  before,  there  is  a  limit  on  the  voltage  that  can  be  applied  to 
the  piezo  material  as  an  actuator.  This  limit  is  due  to  two  factors:  first,  the 
breakdown  voltage  of  the  piezoceramic  and  second,  the  limit  imposed  by  the  sat¬ 
uration  voltage  of  the  power  amplifier.  So,  there  are  constraints  of  the  following 
form  on  inputs, 


Vmin  5:  Vi  2:  Vmax  *  —  1;  *'*  5  Aj  (4.2.5) 

where,  Vi  is  the  input  voltage  to  the  power  amplifier. 

The  Stabilizing  Feedback  Gain  ’L’ 

As  mentioned  in  Step  2,  the  objective  is  to  find  the  feedback  gain  matrix 
L  =  [h  \  h\  ■  ■  ■  | /TO]  such  that  the  input  V  =  LTx  stabilizes  the  system  (asymp¬ 
totically)  while  the  inputs  Vi  satisfy  inequalities  (4.2.5)  for  all  x  €  D  (this 
guarantees  that  DcE  (see  [13])).  Given  the  feedback  gain  L  the  maximum 
of  if x  is  achieved  when  x  is  computed  from  the  following  rule, 
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(4.2.6) 


*£j  max 


I  Mj  Iji  >  0 

I  mj  <  0 

j  =  1,  •••  ,  2 TV. 


Similarly,  the  minimum  of  ifx  is  achieved  for 


% j  min 


=  < 


rrij  Iji  >  0 


Mj  Iji  <  0 

;  =  1,  •••  ,  2AT. 


(4.2.7) 


One  way  to  choose  the  feedback  gain  L  is  to  use  the  LQR  method  and 
increase  the  penalty  on  ipput  until  inequalities  (4.2.5)  are  satisfied  with  Vi 
replaced  by  Ifx ,  for  all  x  €  D.  To  check  if  the  given  feedback  gain  is  acceptable, 
the  above  rules  for  Xj  m,n  and  Xj  max  may  be  used. 


Feasibility  Check 

In  this  step  we  check  to  see  if  the  following  inequality  holds  for  a  given  choice 
of  P  ,  L  ,  and  D  [13]: 

sup  xTPx  <  minxTPx.  (4.2.8) 

x&D  X^SE 

At  the  left  hand  side  of  the  above  inequality,  the  supremum  can  be  achieved 
only  at  the  vertices  of  the  polytope  (a  convex  hull  of  finitely  many  points)  D  . 
The  reason  for  this  is  that  P  can  be  diagonalized  with  an  orthogonal  transfor¬ 
mation  T.  Then  in  the  new  coordinate  system,  the  left  hand  side  of  (4.2.8) 
can  be  rewritten  as,  supzgT(c)  XT'  Z1  (^»  >  0).  The  result  follows  by  noting 
that  first,  the  transformation  induced  by  T  is  actually  a  rotation  and  vertices 
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are  transformed  to  vertices  and  second,  the  maximum  of  the  above  expression 
occurs  at  a  vertex  of  the  transformed  poly  tope.  Therefore,  supx6i;)  xT  Px  can 
be  computed  after  22N  checks  (i.e.  evaluating  xT Px  on  the  vertices  of  the 
region  D  and  taking  the  maximum  value). 

Also,  the  left  hand  side  of  the  inequality  (4.2.8)  is  equal  to  [13], 


rp 

min  x  Px 

x£SE 


min[min{- 


V2- 


V? 


Max 


rfr}]- 


(4.2.9) 


IJP-Hi'  IfP-'U 

Again,  if  the  inequality  (4.2.8)  is  not  satisfied  we  have  two  options:  first, 
change  Q  to  get  another  P  and  second,  increase  the  input  penalty  to  obtain 
another  feedback  gain  L  and  proceed  as  before. 

Gutman  and  Hagander  showed  that  [13]  the  input, 


V  =  sat[(LT  —  kBT  P)x],  (4.2.10) 

stabilizes  the  system  for  all  positive  definite  diagonal  matrices  k.  In  order  to 
choose  the  tuning  matrix  k  we  plot  the  eigenvalues  of  the  closed  loop  system 
(without  the  saturation)  as  k  changes  (here  we  choose  k  —  cl  and  we  vary 
c  to  obtain  the  root  locus).  The  tuning  parameter  k  has  been  chosen  so  that 
the  dominant  eigenvalue  (the  one  that  is  the  closest  to  the  imaginary  axis)  be 
as  far  as  possible  from  the  imaginary  axis.  Figure  4.6  shows  a  portion  of  the 
root  locus. 

The  results  of  simulation  for  this  particular  choice  of  tuning  parameters  are  . 
shown  in  figure  4.7. 

In  order  to  compare  the  performance  of  the  LQR  method  and  the  Gutman- 
Hagander  method,  a  disturbance  in  the  form  of  a  sinusoid  is  applied  at  the  first 
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Root  locus:  K=0.01  -  0.089 


Figure  4.6:  Location  of  the  eigenvalues  of  the  closed  loop  system  as  the  tuning 
parameter  k  varies 

joint.  In  figure  4.8  the  frequency  of  this  sinusoid  is  chosen  to  be  the  lowest 
natural  mode  of  the  beam.  Figure  4.9  shows  the  results  when  the  frequency  of 
the  disturbance  is  the  highest  natural  mode  of  the  beam. 

Up  to  this  point  we  studied  the  behavior  of  the  linearized  system.  To  observe 
the  effect  of  non-linearities(which  come  into  the  picture  when  we  use  the  non¬ 
linear  equations  of  motion,  (2.1.16),  instead  of  the  linearized  ones).  In  order  to 
perform  the  simulation,  the  non-linear  equations  of  motion  has  been  transformed 
in  the  following  form  : 


Qa  =  /(a,  a,  M), 
where  Q  is  a  symmetric  matrix 


(4.2.11) 
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Figure  4.7:  Controlled  vibration  of  the  beam  (Gutman-Hagander  method) 

Two  simulations  has  been  performed.  In  the  first  simulation  we  want  to  show 
the  effect  of  non-linearities  in  the  free  oscillation  of  a  four  link  model.  One  may 
compare  figure  2.5  and  figure  4.10  . 

In  the  second  simulation  the  Gutman-Hagander  control  is  applied  to  the 
linear  and  nonlinear  systems.  Although  in  this  simulation  we  consider  large  de¬ 
flections  of  the  beam  (which  is  unrealistic),  again  comparison  of  figure  4.11  and 
4.12  shows  that  the  responses  are  very  similar  and  the  controller  work  well  even 
for  the  nonlinear  system. 
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Figure  4.8:  Performance  of  Gutman-Hagander  controller  under  external  distur¬ 
bance  (low  frequency  disturbance) 
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Figure  4.9:  Performance  of  Gutman-Hagander  controller  under  external  distur¬ 
bance  (high  frequency  disturbance) 
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Figure  4.10:  Free  Oscillation  (Nonlinear  Model) 
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Figure  4.11:  Controlled  vibration  of  the  beam  (Large  deflection,  Linear  Model, 
G-H  method) 
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Non-Linear  Model 
k=0.022I 


Figure  4.12:  Controlled  vibration  of  the  beam  (Large  deflection,  Nonlinear 
Model,  G-H  method) 
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Chapter  5 


Observer  Based  Gutman-Hagander 
Type  Controller 


In  this  Chapter  we  prove  that  under  certain  conditions,  a  Gutman-Hagander 

* 

type  controller  with  a  classical  state  observer  can  be  used  for  LTI  systems  with 
constraints  on  the  magnitude  of  input  to  make  them  stable. 


5.1  Introduction 

The  problem  of  stabilizing  LTI  (linear  time  invariant)  systems  with  constraints 
on  the  magnitude  of  the  input  has  been  considered  by  Gutman-Hagander  [13]. 
The  method  requires  the  knowledge  of  all  the  state  variables.  In  many  practi¬ 
cal  situations  it  is  not  possible  or  it  is  too  expensive  to  measure  all  the  state 
variables.  So,  one  may  try  to  use  the  estimated  values  of  the  state  variables 
and  then  exploit  the  Gutman-Hagander  (G-H)  type  controller.  We  prove  that 
it  is  possible  provided  that  some  conditions  are  satisfied.  In  Section  5.2  we 
obtain  a  proper  Lyapunov  function  to  prove  the  stability  of  the  method.  In 
Section  5.3  we  propose  an  algorithm  (which  is  very  similar  to  the  one  proposed 
by  Gutman  and  Hagander)  to  choose  the  control  parameters.  Next,  in  order  to 
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test  the  performance  of  the  control  strategy  we  apply  this  control  law  to  the 
problem  of  damping  the  vibration  of  a  cantilever  beam. 


5.2  Control  Strategy 


Consider  a  controllable  and  observable  linear  time  invariant  plant  represented 
as: 


x 


< 


y 


=  Ax  +  Bu, 
=  Cx. 


(5.2.1) 


The  control  strategy  due  to  Gutman  and  Hagander  proposed  the  following 
control  law  as  stabilizing  input: 


u  =  sat[{LT  -  I<BtP)x\, 


(5.2.2) 


where  if  is  a  positive  definite  diagonal  matrix.  In  what  follows  we  prove  that  the 
above  control  law  still  makes  the  system  stable,  even  if  one  uses  the  estimated 
state  variable,  x  ,  computed  by  an  state  observer  instead  of  the  actual  sate,  x. 
i.e. 


u  =  sat[{LT  -  I< BT P)x\.  (5.2.3) 

Let  e  be  the  estimation  error  (i.e.  e  =  x  —  x).  It  is  well  known  that  the 
dynamics  of  the  error  is  described  as: 

e  —  (A  —  HC)e,  (5.2.4)' 

where,  H  is  the  observer  gain. 
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Let  P  ,  G  be  the  positive  definite  solutions  to  the  algebraic  Lyapunov 
equations  : 

PAC  +  A?P  =  -Q,  (5.2.5) 

GA0b  +  AjbG  —  —M.  (5.2.6) 

Where  Aob  =  A  —  DC ,  Ac  =  A  +  LT B  ,  and  M  and  Q  are  symmetric 
positive  definite  matrices  M  (these  positive  definite  solutions  exist  because 
Ac  and  A0f,  are  stable).  Suppose  R  is  a  positive  definite  diagonal  matrix  and 
apply  the  following  input  to  the  plant  5.2.1: 

u  =  (Lt  —  RBT  P)x  (5.2.7) 

The  claim  is  that  V^x,t)  defined  below  is  a  Lyapunov  function  for  the 
closed  loop  system  (including  the  observer). 


V  =  xTP  x  +  aeTGe  (5.2.8) 

for  some  positive  constant  a  which  is  to  be  chosen  to  make  the  derivative  of 
V  negative.  One  can  show  [13]  that  , 

sat[(LT  -  KBtP)x]  =  ( Lt-RBtP)x ,  (5.2.9) 

for  0  <  R  <  K,  provided  that  the  condition  D  C  E  holds  (see  section  4.2 
for  definitions  of  D  and  E ).  Taking  the  derivative  with  respect  to  time  of 
eq.  5.2.8  and  substituting  for  u(t)  from  eq.  5.2.7  we  get: 
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V  =  (Ax  +  Bu)TPx  +  xtP(Ax  -f  Bu)  +  aeTA?bGe  +  aeTGA0i,e 


—  xT(PAc  +  A^P)x  —  xT  PRBT  Px  —  xT  PBLTe 
+xTPBRBTPe  -  eTLBTPx  -  xTPBRPx 
+eT  PB  RBT  Px  +  aeT  A^bGe  +  aeTGA0be 


(5.2.10) 


=  —xTQx  —  2xT(BTP)TR(BTP)x 
-2 er(LBTP  -  PBRBTP)x  -  aeTMe. 

To  simplify  the  above  equation  define  S  and  U  as, 


S  =  QA2(BTP)TR(BTP),  S  >  0, 
U  =  2LBtP  -  2PBRBtP. 

Rewrite  5.2.10  in  terms  of  S  and  U  to  get: 


(5.2.11) 


V  =  -xTSx-eTUx-aerMe.  (5.2.12) 

The  idea  is  to  choose  a  large  enough  such  that  V  becomes  negative.  Since 
matrices  M  and  S  are  symmetric  and  positive  definite  these  inequalities  hold: 


\™n\\e\\2<eTMe<X%?*\\e\\\ 

(5.2.13) 

Ar||e||2  <  eTSe  <  A£aI||e||2, 

(5.2.14) 

eTUx<\eTUx\<\\U\\\\x\\\\e\\. 

(5.2.15) 

Thus, 
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v  <  -Aji"w:'-oAsrii«r+ ni/ii wwi 

=  -2a|||x|P  -  (g)2||e||2  +  ||C/||||e||l|x||  +  (g)2||e||!  -  aAjf  Nit 

where,  2e*|  = 

Subsequently, 


V  <  -a|||x||2.  (5.2.18) 

5.3  Stability  Under  Variable  Matrix  R 

In  the  previous  section  we  proved  the  stability  of  the  plant  when  the  control 
5.2.7  is  used.  To  relate  this  control  to  the  control  with  saturation  given  by 
Equation  (5.2.3)  (see  the  inequality  (5.2.9)),  we  need  the  following  proposition. 
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Proposition:  The  input  5.2.7  with  the  condition  0  <  R  <  K  makes  the 
plant  5.2.1  asymptotically  stable. 

Proof: 

Now,  both  A™1"  and  \\U\\  in  5.2.17  are  variables.  So,  we  try  to  obtain  an 
upper  bound  for  \\U\\  and  a  lower  bound  for  A™”.  From  5.2.11  we  can  write: 


\\U\\  =  \\2LBtP  -2PBRBtP\\ 

<  \\2LBtP\\+2\\PB\\2\\R\\  (5.3.1) 


<  l|2IBTP||  +  2||PB|n|/C||, 

* 

the  last  inequality  follows  from  the  assumption  that  R  <  K. 
The  next  step  would  be  to  find  an  upper  bound  for  Xgm: 

A™'*"  =  inf||3.||=1  xTSx 

=  infj|a;||-1(a:TQx  +  xTNx). 
where  N  =  2 (BT P)T R(BT P).  We  also  have: 

Ag’71  <  xTQx 

<  xTQx  +  xTNx. 


(5.3.2) 


(5.3.3) 


The  above  inequality  shows  that  A g1"  is  a  lower  bound  for  the  set  defined 
by  {xT Qx  +  xTNx  |  ||x||  =  l|  we  get: 


ymm  ^  ymtn 


(5.3.4) 
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Therefore,  the  inequality  5.2.17  is  satisfied  if  the  following  holds: 


a  > 


2\\LB’rP\\+2\\PBni<\\ 


\min 

aq 


(5.3.5) 


This  implies  the  negativity  of  V  and  asymptotic  stability  of  the  plant  repre¬ 
sented  by  5.2.1  with  the  control  given  by  Equation  5.2.3. 


5.4  Controller  Design  Procedure 

The  controller  design  algorithm  is  essentially  the  same  as  the  one  proposed  by 
Gutman-Hagander.  The  only  difference  is  that,  now  we  have  to  redefine  the  sets 
D  ,  E  ,  etc.  introduced  in  [13].  Lets  define  the  augmented  state  variable,  x  ,and 
the  matrix  P  as, 


x 


(5.4.1) 


P  = 


P  0 
0  olG 


(5.4.2) 


Now,  V  can  be  written  in  terms  of  the  new  state  variable: 


V  =  xTPx.  (5.4.3) 

To  obtain  the  various  control  parameters  we  follow  these  steps: 

Step  1:  Find  the  set  of  initial  conditions,  D  ,  as  restrictive  as  possible. 
This  set  must  be  large  enough  to  include  all  the  possible  initial  values  of  the 
augmented  state  variable  x  that  may  occur  in  practice. 
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Step  2:  Choose  the  stabilizing  feedback  gain  L  in  a  way  that: 

g  <  17 x  <  h  \/x  €  jD,  (5.4.4) 

where,  L  =  [I  —  I]T L  and  L  is  a  stabilizing  feedback  gain  for  the  plant 

5.2.1  and  h,  g  €  Rm  are  the  upper  bound  and  the  lower  bound  on  the  input 
respectively. 

Step  3:  Choose  the  diagonal  matrix  K  with  positive  elements  and  solve  the 
Lyapunov  equation  5.2.5  for  some  positive  definite  matrix  Q.  Then  choose  a 
positive  constant  a  according  to  the  inequality  5.3.5. 

Step  4:  Choose  some  positive  definite  matrix  M  and  solve  the  Lyapunov 
equation  5.2.6  to  get  the  matrix  G.  Construct  the  matrix  P  according  to 

5.4.2  and  change  the  control  parameters  until  the  following  inequality  holds: 

sup  xT  Px  <  mina:TPx.  (5.4.5) 

xeD 

Step  5:  Compute  the  estimated  state,  x,  using  a  state  observer  and  apply 
the  stabilizing  input  5.2.7  to  the  plant. 

5.5  Simulation 

To  justify  the  second  step,  note  that  by  the  separation  theorem,  if  u  =  17 x 
stabilizes  the  plant  5.2.1,  the  input  u  =  LTx  stabilizes  the  whole  system  in¬ 
cluding  the  observer  (provided  that  Aab  is  stable  too).  Thus,  the  procedure  for 
computing  an  acceptable  LT  is  to  find  stabilizing  L  ,  say  by  LQ  method,  and 
then  define  L  as  above.  Note  that  we  can  write: 

lTx  =  Tt£,  (5.5.1) 
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this  means  that: 


x  G  E  <=>  x  €  E,  (5.5.2) 

with  the  obvious  definitions  for  E  and  E. 

The  implication  5.5.2  is  important  because  in  the  proof  of  the  stability  we 
rely  on  the  fact  that  x  remains  in  the  set  E  and  by  5.5.2  we  know  that  this 
is  guaranteed  provided  that  the  following  holds: 

D  C  n  C  E.  (5.5.3) 

Again,  we  can  increase  Jihe  input  penalty  so  that  L  satisfy  5.4.4.  To  see  how 
efficient  the  proposed  control  strategy  is,  we  apply  the  method  to  the  problem 
of  damping  the  vibration  of  a  cantilever  beam.  We  assume  that  the  deflection 
angles  are  directly  measurable  and  we  design  an  observer  to  construct  the  other 
state  variables  (time  derivative  of  the  deflection  angles). 

5.5.1  Observability 

Consider  the  state  variable,  x  defined  as  the  vector  of  alternating  deflection 
angles  and  their  respective  corresponding  time  derivative: 
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£*1 

di 


X 


Cti 

di 


Oin 

d-n 


(5.5.4) 


It  can  be  shown  that  the  A  and  C  matrices  describing  the  cantilever  beam 
are  of  the  form: 


A  = 


0  1 
*  0 
0  0 
*  0 


0  0 
*  0 
0  1 
*  o 


0  0 
*  0 
0  0 
*  0 


0  0  0  0  0  1 

*  0  *  0  *  0 


(5.5.5) 


1  0  0  0  0  0 
0  0  1  0  0  0 


0  0  0  0  1  0 


(5.5.6) 


here,  for  the  sake  of  simplicity,  we  showed  A  and  C  for  n  =  3  (  *  stands  for 
nonzero  elements). 

Therefore,  the  observability  matrix  has  the  following  form: 
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o  = 


1  0  0  0  0  0 

0  0  1  0  0  0 

0  0  0  0  1  0 

0  1  0  0  0  0  , 

0  0  0  1  0  0 

0  0  0  0  0  1 


t 


(5.5.7) 


which  is  full  rank  and  consequently  the  pair  (C,A)  is  observable. 


5.5.2  Simulation  Results 

The  observer  settling  time  is  chosen  to  be  much  faster  than  that  of  the  closed  loop 
system.  Note  that  the  faster  you  make  the  observer  the  larger  the  estimation 
error  would  be  and  the  latter  cause  a  larger  initial  condition  set  D.  Thus  in 
a  sense  there  is  a  trade  off  between  the  size  of  the  set  D  and  the  speed  of  the 
observer.  Figure  5.1  shows  the  estimation  error  of  the  observer  for  a  particular 
choice  of  the  parameters: 

Now,  we  can  set  up  the  controller-observer  according  to  the  scheme  shown 
in  figure  5.2. 

Although  the  analysis  has  been  performed  for  the  linear  systems,  simulation 
results  show  that  the  response  of  the  linear  and  nonlinear  systems  are  very 
similar. 
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Figure  5.1:  Estimation  error  in  deflection  and  velocity 

Shown  in  figure  5.3  is  the  response  of  the  nonlinear  model  for  the  cantelevered 
beam  when  the  scheme  shown  in  figure  5.2  has  been  used  to  damp  the  vibrations. 


Figure  5.2:  G-H  type  controller  with  observer 
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Input  voltage  [volts] 


Time  [sec] 
Saturated  Input 


Time  [sec] 


Non-Linear  Model 

k=0.0I2I 


Figure  5.3:  G-H  type  controller  with  observer 


Chapter  6 


Numerical  Solution  to  the  PDE  for  the 
Cantilever  Beam  With  Control  and 
Derivation  of  the  Closed  Form 
Transfer  Function 


The  purpose  of  this  chapter  is  to  first  develop  a  simulation  tool  to  study  the 
behavior  of  the  PDE  describing  the  transversal  deflection  of  an  Euler-Bernoulli 
beam  and  then  to  exploit  it  to  derive  an  empirical  transfer  function  matrix  of  a 
set  up  that  is  used  to  damp  out  the  vibrations  of  the  beam.  To  solve  the  PDE 
(numerically)  an  implicit  finite  difference  method  is  used.  The  approximation 
is  done  using  a  scheme  based  on  least  squares  proposed  in  [15],  that  employs 
rational  wavelets.  A  closed  form  expression  for  the  transfer  function  of  the 
system  is  obtained  at  the  end  and  comparison  has  been  made  with  the  empirical 
transfer  function  estimated  by  the  spectral  analysis. 
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6.1  Introduction 


The  transversal  deflection  of  a  beam  can  be  described  by  a  partial  differential 
equation  and  proper  boundary  conditions  which  reflect  the  type  of  constraints 
at  its  ends. 

In  the  process  of  designing  a  controller,  to  damp  out  the  vibrations  of  the 

beam,  one  may  fit  a  finite  dimensional  model  to  the  beam  or  use  some  distributed 

parameter  controller.  In  either  case,  to  test  the  efficiency  of  the  controller  we 

need  to  apply  the  control  method  to  a  model  of  the  beam,  which  takes  into 

account  the  infinite  dimensional  nature  of  the  system.  Here,  we  developed  a 

simulation  program  to  solve  the  PDE  pertaining  to  Euler-Bernoulli  beam.  At 

first,  we  had  considered  the  free  oscillation  of  the  beam,  and  then  to  verify  the 

* 

results,  the  total  energy  of  the  beam  as  a  function  of  time  was  computed  (since 
no  damping  mechanism  exists  this  should  be  a  constant).  Next,  we  include  the 
effect  of  internal  and  viscous  damping  in  the  model  to  make  it  more  realistic. 
With  the  simulation  tool  in  hand  we  apply  a  swept  sine  wave  as  the  input  to  the 
actuating  piezo  material  and  observe  the  voltage  across  the  piezoceramics  that 
are  used  as  sensors.  The  input  and  output  waveforms  then  enable  us  to  get  an 
empirical  frequency  responce.  In  the  next  step  we  fit  a  finite  dimensional  model 
to  the  frequency  response. 


6.2  Finite  Difference  Approximation 

The  following  PDE  describes  the  transversal  vibrations  of  a  thin  beam  [1], 


d2y  ,  dy  d2  ,T71Td2y 


dt 2 


+  -<m  + 


dPy 


0,  (6.2.1) 
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where,  y  =  y{x\  t)  is  the  displacement  from  the  equilibrium  position  at  distance 
x  along  the  beam  from  one  end,  7  is  the  coefficient  of  viscous  damping,  p  is  the 
linear  mass  density,  and  C£>I  is  a  coefficient  reflecting  the  internal  damping. 
Define  two  new  variables,  v  and  w  as  below 


Vt 


w  = 


lyxx 


(6.2.2) 


XX  t 


Now,  PDE  (6.2.1)  can  be  written  in  the  following  form 


vt 


'-v  —  awxx 


Wt  —  avxx  -pjfg jVt  +  -jjgp[Vxxt 


(6.2.3) 


where  a  =  \J El / p  and  subscripts  x  and  t  are  stand  for  partial  derivatives  with 
respect  to  x  and  t. 

One  possible  finite  difference  representation  of  the  above  PDE  is 


<+1~< 

At 

w”+1-w 

Zl _ l 


(SW+1+(S2wn 
~  a  2(Az)2 


>1 


1  -  (a  _  _ cp  y  1  zcpj 

~  v°  jmD  UAxY  Vt  + 


(62v)2+1-(62v)? 


(6.2.4) 


At  ~  V-  2(A^P  7pElrt^VpEl  (Ar)2(A  t) 

where,  superscripts  and  subscripts  represent  time  and  position  indices,  respec¬ 
tively  and  is  the  finite  difference  approximation  of  the  second  derivative  of 
w  with  respect  to  x.  That  is, 


(62w)?  =  wj+1-2w]  +  wj_1.  (6.2.5) 

Of  course,  instead  of  (6.2.4),  one  could  use  some  other  scheme  of  finite  dif¬ 
ference  approximation,  but  the  advantage  of  using  (6.2.4)  is  that  it  makes  the 
difference  equation  unconditionally  stable  (i.e.  stability  does  not  depend  on  the 
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choice  of  step  sizes  Ax,  and  At).  Nevertheless,  system  (6.2.4)  entails  other  com¬ 
plications  due  to  its  implicit  nature,  and  the  extra  effort  required  to  solve  the 
system  (6.2.4). 

We  may  cast  (6.2.4)  into  the  form 

-Auj+1  +  Buj  -  Auj_x  =  dj,  j  =  1,  2,  •  •  •  ,  J  -  1  (6.2.6) 

in  which, 


a  — 


aAt 
(Ax)2  * 


P  =  2[( 


-  f  1  — r~JTT . !  cnla _ 

s/pEI '  2(Aa;)2  '  \/pEI(Ax)2 


cp4f_T/l 


Ui  = 


n+1 


W 


n+1 


A  = 


£  o 
2  U 


(6.2.7) 


B  = 


1  ~a 

P  1 


dj 


-  !K+1  “  2 w]  +  «£_,)  -  ^ 
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6.3  Preliminaries 


In  the  algorithm  that  we  are  going  to  use,  it  is  required  that  the  boundary 
conditions  to  be  specified  in  a  form  shown  below  (in  the  next  section  we  show 
how  to  write  down  these  equation  for  the  cantilever  beam) 

uji  =  Huj1+1  +  l,  (6.3.1) 

uj—j^—x  =  Muj-j2  +  n,  (6.3.2) 

with  known  matrices  M ,  H  and  vectors  n  and  1. 

Now,  consider  the  difference  equation, 


Uj  —  EjUj+i  +  fj,  j  —  0,  1,  •••  (6.3.3) 

* 

Equations  (6.2.6)  and  (6.3.2)  specify  a  two-parameter  family  of  solutions  of 
(6.2.6).  Equation  (6.3.3)  also  has  a  two-parameter  family  of  solutions  and  these 
two  families  are  identical  if  we  choose  E  and  /  as, 

Eh  =  H, 

fh  =  l'  (6.3.4) 

Ej  =  (B  -  AEj.^-'A, 

fj  =  (B  —  AEj-i)~1(dj  +  Afj-i)  j  >  1. 

Rewrite  (6.3.3)  for  the  index  J  —  jh  —  1  to  get 

uj,_j2-i  —  Ej-j2-xUj-j2  +  fj—j2-\.  (6.3.5) 

Eliminating  uj-j2- 1  between  (6.3.5)  and  (6.3.2)  give  us  uj-j2  in  terms  of  . 
other  variables: 


uJ-h  =  {Ej-j2-i  —  M)  1(n  —  /j_j2_i).  (6.3.6) 
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Now  that  we  have  defined  the  various  variables  involved,  we  can  state  the 


algorithm. 


6.4  Algorithm 

Assume  that  the  state  of  the  system,  u”,  j  =  0,  1,  •  •  •  ,  J,  for  some  n  is  known. 
The  algorithm  described  below  shows  how  to  obtain  u”+1,  j  —  0,  1,  •  ■  •  ,  J 
at  the  next  time  step. 

•  Step  1:  Solve  the  difference  equations  (6.3.4)  for  j  =  j  i,  ji  +  1,  •••  ,J  —  1 
with  the  initial  values  specified.  (For  the  problem  that  we  want  to  study, 
it  can  be  shown  that  Ei  and  /,■  both  remain  bounded) 

4. 

•  Step  2:  Using  the  matrix  and  the  vector  1  compute  uj_;-2 

from  equation  (6.3.6). 

•  Step  3:  Compute  u,-,  i  —  J  —  j?  —  1,  •  •  •  ,  ji  +  1,  j i  by  solving  (6.3.3) 
backward  in  j. 

•  Step  4:  ft  remains  to  obtain  a  few  of  the  state  variables,  u,-,  at  each  end. 
i.e.  u0,  •  •  •  ,  Ujj- 1  at  one  end  and  uj,  •  •  • ,  uj_J2+1  at  the  other  end.  These 
can  be  obtained  from  (6.2.6) 

6.5  Boundary  Conditions  for  the  Cantilever 
Beam 

In  order  to  use  the  above  algorithm  we  need  to  specify  the  boundary  conditions 
for  the  cantilever  beam  in  the  form  of  equations  (6.3.2)  and  (6.3.2).  So,  the 
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parameters  ji,  j2)  M,  H,  n ,  and  l  are  need  to  be  computed.  The  boundary 
conditions  are: 


y(t,  0)  =  yx(t,  0)  =  0  Vi  >  0  (clamped  end) 

M(i,L )  =  Mx(t,L)  =  0  Vi  >  0  (free  end) 

where,  M  is  the  total  moment  (internal  and  external)  defines  by: 


(6.5.1) 


+  (6.5.2) 

The  boundary  condition  at  the  clamped  end  implies  that  Vo  =  iq  =  0. 
Equation  (6.2.6)  together  with  the  boundary  conditions  at  the  clamped  end  have 
a  solution  which  is  linearly  dependent  on  two  parameters.  Choose  A  =  [wq  uq]T 
as  the  vector  of  parameters.  Then  one  can  write 


uj  =  UjX+gj,  j  =  0,  1,  •  •  •  ,  J,  (6.5.3) 

where,  Uj  is  defined  by, 

OTU  =  |f-  (6-5.4) 

The  procedure  is  to  compute  Uj  for  i  =  0,  1,  •  •  •  ,  j\  such  that  Uj1  is  the 
fist  invertible  matrix  in  that  sequence  of  matrices.  Then  we  have 


^  —  Uji  {uji 


(6.5.5) 


which  in  turn  results, 


uii— i  —  Ujl—\Ujl  (iij1  gjf)  +  gj1  -i 


(6.5.6) 


82 


which  is  in  the  form  of  (6.3.2)  with  H  and  l  as  below, 


H  =  Uh 


(6.5.7) 


i  —  +  9h-i- 


For  the  problem  considered  here  f/2  is  the  first  nonsingular  one  and  it  can 
be  readily  verified  that  the  previous  ones  are, 


U0  = 


0  1 
1  0 


(6.5.8) 


Ui  = 


0  0 
0  1 


(6.5.9) 


U2  = 


0  1 
u  0 


-1  2 

Now,  we  need  to  compute  gi  and  g2: 


(6.5.10) 


<72  —  U2  —  U2\ 

=  A~lBui  —  u0  —  A_1di  —  U2X  (6.5.11) 

Also,  for  <71  we  have, 


gi  =  «i  -  U\X 
0 
0 


Thus, 


(6.5.12) 


A  =  2,  H  =  IhU*1,  l  =  IhUZ'A-'dLi.  (6.5.13) 
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Now,  consider  the  boundary  conditions  at  the  free  end  of  the  beam.  Again, 
from  (6.5.1)  we  get  wj  ---  wj-\  =  0  and  this  time  we  choose  A  =  [vj  vj-\]T  as 
the  parameters.  Similarly,  define  Uj  by, 


(j  n  _  ditji 
Kuj)ik 


(6.5.14) 


8vj-k+ ! 

Going  backward  in  position  from  the  free  end  toward  the  clamped  end,  Uj- 2 
is  the  first  nonsingular  matrix  and  we  have, 


Uj 


1  0 
0  0 


(6.5.15) 


Uj- 1  = 


0  1 
0  0 


(6.5.16) 


Uj- 2  = 


-1  2 
0  =* 


(6.5.17) 


Thus,  j 2  =  2.  Following  the  same  steps  as  for  the  clamped  end  but  this  time 
for  the  free  end,  one  gets, 


^  =  Uj_2(uj-2  —  gj-2), 


(6.5.18) 


and, 


uj~  3  —  U j —3U j  ^_2{.u  j —2  —  gj- 2)  +  gj- 3- 


(6.5.19) 


Again,  it  can  be  shown  that  gj-2  =  —A  1dj_1.  To  obtain  gj- 3  first  we 
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need  to  compute  t/j_3.  So,  write 


uj- 3  =  A  lBuj- 2  -  uj-i  -  A  1dj_2, 

—  A~l B(A~1Buj- 1  —  uj  —  A-1dj_i)  -  uj_ i  -  A_1dj_2, 


q _ 4_  SS 

c*/3  /3 


a/3 


uj_!  —  A  1Buj  —  A  XBA  1dj_1  — 


(6.5.20) 


A  1dj_2. 


from  which  we  can  easily  obtain  Uj- 3: 


t/j- 


-2  3-^ 

1  ___  8. 

a  a 


(6.5.21) 


Substitute  for  uj_3  and  Uj-z  into  the  relation  for  #j_3  and  simplify  to  get 


Qj- 3  =  nj-3  —  Uj- 3  A, 

(6.5.22) 

=  -A^BA^dj-i  -  A~xdj-2- 

To  compute  M  and  n,  plug  in  the  expressions  that  we  obtained  for  gj-3, 
gj- 2,  U j—zi  and  Uj-2  into  (6.5.19) 

M  = 

n  =  (6-5.23) 

A-1BA-ldj-1-A-1dJ_2. 

So,  we  have  computed  all  the  terms  that  are  required  by  the  algorithm. 


6.6  Simulation 

To  implement  the  algorithm  mentioned  above,  first  we  must  specify  the  initial 
shape  and  velocity  of  the  beam.  We  assumed  a  uniform  load  on  the  beam  of 
length  1  m  and  the  load  is  such  that  it  results  in  a  deflection  of  8cm  at  the  tip. 
The  deflection,  y,  at  t  —  0  and  at  distance  x  from  the  clamped  end  is  [5], 
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y  =  -^(QL2x2 -4Lx3  +  x4).  (6.6.1) 

where,  8  is  the  maximum  deflection  (which  occurs  at  the  tip)  and  L  is  the  length 
of  the  beam. 

First,  we  set  both  the  internal  and  damping  coefficients  to  zero.  This  gives 
us  a  way  to  verify  the  correctness  of  the  algorithm.  Because  in  the  absence  of 
damping,  the  total  energy  must  remain  constant  with  time.  Figure  (6.1)  shows 
the  velocity  and  deflection  of  various  fixed  point  on  the  beam  versus  time. 


Figure  6.1:  Velocity  and  deflection  as  a  function  of  time  (no  damping) 


The  velocity  and  deflection  at  the  free  end  of  the  beam  are  shown  in  figure 

(6.2) 
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Figure  6.2:  Velocity  and  deflection  of  the  free  end  as  a  function  of  time 

Also,  we  can  show  the  shape  of  the  beam  as  it  oscillates.  Each  curve  in  figure 
(6.3)  represents  the  shape  of  the  beam  at  a  fixed  instant  of  time. 

The  energy  of  the  system  at  time  t  is  obtained  from  the  following  relation: 


<6-6-2> 

Using  the  first  order  approximation  of  the  integral  the  above  equation  can 
be  easily  written  in  terms  of  the  variables  v  and  w  introduced  at  the  beginning 
of  the  chapter, 


En 


EMW)2)  +  K)2))- 


i= 0 


(6.6.3) 
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Shape  of  the  beam  at  different  instants  of  time 


Figure.  6.3:  Free  oscillation  of  the  beam 

Since  the  deviation  of  the  energy  from  its  average  value  is  very  small,  in 
figure  (6.4)  we  show  the  percentage  of  deviation  of  it  versus  time. 

As  it  can  be  seen,  the  error  is  relatively  small  and  perhaps  the  major  part 
of  the  error  is  due  to  the  first  order  approximation  that  we  used  for  integration 
of  (6.6.2).  Introducing  internal  and  viscous  damping  give  us  the  following  set  of 
graphs. 


6.7  Estimating  the  Empirical  Transfer  Func¬ 
tion 

We  are  interested  in  finding  the  transfer  function  from  the  output  which  is  the 
voltage  across  piezoceramics  attached  to  the  beam  and  operate  as  sensors,  to 
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Percentage  of  the  variation  of  the  total  energy  vs  time 


Figure  6.4:  The  energy  stored  in  the  beam  remains  constant. 

the  input  voltage  that  is  applied  to  the  actuating  piezoceramics.  The  model 
that  we  use  is  shown  in  Fig.  (6.7). 

Since  the  PDE  (6.2.1)  is  linear,  we  may  use  the  ratio  of  Fourier  transform  of 
the  output  to  that  of  the  input  and  this  give  an  estimate  of  the  transfer  function 
[9].  That  is, 

<9(e”)  =  —El  (6.7.1) 

UN(u>) 

and  Yn(u)  is  defined  by, 

VVH  =  (6.7.2) 

It  can  be  shown  that  [9]  the  following  relationship  holds  between  the  empir¬ 
ical  transfer  function,  (^(e-7"),  and  the  true  transfer  function  G(eJUJ): 

G(e^)  =  G(e?w)  +  77~7~t>  (6-7.3) 

Ui<r{LO) 
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8  10  12 
Time  [sec] 

Deflection  of  the  tip  vs  time 


Time  [sec] 

Figure  6.5:  Velocity  and  deflection  as  a  function  of  time  (with  passive  damping) 
where  Rn(w)  decays  as  1/y/N. 

In  the  case  of  periodic  input  the  above  procedure  gives  a  good  estimation  of 
the  transfer  function  around  the  frequencies  that  exist  in  the  periodic  input.  But, 
if  the  input  is  not  periodic  we  must  assume  that  the  values  of  the  true  transfer 
function  are  related  at  different  frequencies.  As  an  alternative  procedure,  we 


may  obtain  the  spectral  estimate  GTeJU')  from  the  expression, 


BE 


x  10'2  Total  energy  vs  time 


Figure  6.6:  The  energy  stored  in  the  beam  decreases  with  time  as  a  result  of 
passive  damping 

where, (u>)  is  the  power  spectrum  of  the  input  and  is  obtained  from, 

OO 

T~  —  OO 

with  R ^  (r)  as  the  autocorrelation  function  defined  as, 


Figure  6.7:  The  inputs  and  outputs  to  the  system 
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Ru(t)  = 

JV  t=l 

and  ro7(r)  is  a  lag  window  with  shape  parameter  7  and  is  chosen  such  that 
vanishes  for  |r|  >  where,  S1  «  N.  This  allows  us  to  approximate  6.7.5 
with  [9]. 

$JV)  =  £  m,(r)^r(r)e-’>". 

T  —  —  6-y 

The  expression  for  is  similar.  There  are  different  choices  for  the  lag 

window  u;7(r)  among  which  we  use  the  Hamming  lag  window  function  defined 
as  [9]. 


1  7 TT 

w^T)  =  -{1  + cos -^)  \t\  <  7,  (6.7.5) 

here,  the  scaling  parameter  7  is  chosen  to  be  equal  to  5^  in  6.7.5.  The  effect  of 
choosing  a  relatively  large  is  that  it  results  in  a  large  bias  in  the  estimated 
transfer  function,  whereas  decreasing  8. y  results  a  larger  variance  in  G(e^). 

To  get  a  “good”  estimate  of  the  transfer  function  the  test  input  should  be 
as  “informative”  as  possible.  A  good  choice  may  be  a  white  noise  because  its 
spectrum  ideally  contains  all  frequencies.  However,  the  input  voltage  to  the 
piezoceramics  ,  H(x,i),  must  be  differentiable  (with  respect  to  time)  as  it  can 
be  seen  from  6.2.1.  In  fact,  piezoceramics  behave  like  small  capacitors  and  a 
discontinuity  of  the  voltage  across  them  theoretically  requires  an  infinite  current. 
A  sweeping  sinusoid  (frequency  modulation  of  a  triangular  waveform)  seems  to 
be  a  reasonable  input  because  it  is  smooth  enough  and  more  importantly  covers 
a  wide  range  of  frequencies. 
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To  demonstrate  the  effectiveness  of  the  method  and  also  to  verify  the  men¬ 
tioned  procedure,  we  used  the  above  method  to  obtain  the  frequency  responce 
of  a  known  system.  Note  that  here  we  are  obtaining  the  frequency  responce  just 
by  one  experiment  and  we  do  not  have  to  wait  until  the  settling  time  of  the 
system.  The  following  system  has  been  used  as  the  test  model, 


H(s) 


289. 

s2  + 0.5s +  289.' 


(6.7.6) 


As  input,  we  excited  the  system  with  the  frequency  modulated  of  a  sinusoid 
that  varies  between  zero  and  uimax.  i.e. 


u(t)  —  sin  ( u(t)t ), 


(6.7.7) 


with, 

oj(t)  =  <+naI(l  —  coscuot)/2.  (6.7.8) 

The  input-output  waveforms  together  with  the  empirical  frequencies  are 
shown  in  figures  (6.8)  and  (6.9). 


The  same  procedure  now  can  be  use  to  form  the  empirical  transfer  function 
of  the  beam.  Here  we  have  two  inputs  and  two  outputs.  Since  the  superposition 
principle  holds  it  is  a  meaningful  thing  to  talk  about  the  matrix  transfer  function. 
In  the  set  up  shown  in  figure  (6.7)  there  are  two  inputs  and  two  outputs. 

The  outputs  are  designed  to  be  the  voltages  that  are  induced  in  the  piezoce¬ 
ramics  due  to  bending  of  the  beam.  Figure  (6.10)  shows  a  portion  of  the  flexible 
beam  which  is  bended.  Assume  that  the  piezo  material  is  attached  to  the  up¬ 
per  (or  lower)  surface  of  the  beam  and  the  objective  is  to  compute  the  change 
of  length  in  it.  The  neutral  axis  is  shown  with  the  dotted  line  and  its  length 
would  not  change  due  to  bending  since  axial  forces  has  not  been  considered. 
The  change  in  slope  for  the  neutral  axis  is, 
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Time[sec] 

Output 


Figure  6.8:*The  inputs  and  outputs  to  the  system 


A6  =  yxx(x)dx. 


With  c  the  distance  from  neutral  axis  to  the  outer  surface  and  from  the  equal 
angles  we  have, 


g 

-  =  A  6  =  yxxdx. 
c 

Thus,  the  strain  in  the  outer  surface  is 


e(x)  = 


e 

dx 


cyxx(z')- 


The  total  change  of  length  in  a  patch  with  its  ends  at  xx  and  x2  is  obtained 


Amplitude 


freq.  [rad/sec] 


Phase 


freq.  [rad/sec] 


Figure  6.9:  Bode  plot  of  the  true  system  together  with  the  estimated  ones 
by  integration, 

A  L  =  f:;e(a)da 

=  Q  cyxx(a)da- 
-  c(yx(x 2)  -  yx(x  1)). 

The  induced  voltage  in  the  piezoceramic  is  proportional  to  the  total  change 
in  length  and  the  constant  of  proportionality  if  called  d^i 

Vout  —  d^\^S.L 

(6.r.9) 

=  d31c(yx(x2)  -  yx(x  1)). 

Shown  in  figure  (6.11)  is  the  Bode  plot  obtained  by  the  spectral  analysis  . 
described  above  (Eq.  6.7.4). 

Since  the  damping  that  has  been  introduced  into  the  model  is  fairly  small,  we 
still  may  expect  that  the  eigenvalues  of  the  system  do  not  change  considerably. 
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Figure  6.10:  The  strain  induced  in  piezo  material 

The  eigenvalues  of  the  undamped  system  can  be  computed  analytically  [8]  and 
the  first  few  of  them  are,  * 


where, 


(6.7.10) 


Ai  =  1.875 /L, 

A2  =  4.695/L, 

A3  =  7.855 /L, 

A4  =  10.996 /L, 

A5  =  14.137/L, 

Based  on  the  above  figures  and  the  parameters  used  in  the  first  fifth  simula¬ 
tion  the  natural  modes  are, 


wx  =  16.68,  cu2  =  104.57,  cj3  =  292.83, 

w4  =  573.852,  w5  =  948.516, 
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mag. 


Figure  6.11:  Bode  plot  of  the  empirical  transfer  function  from  U\  to  y\ 

which  closely  match  the  frequencies  at  which  the  peaks  of  the  modulus  of  the 
transfer  function  in  Fig.  (6.11)  occurs. 

6.8  Fitting  a  Finite  Dimensional  Transfer  Func¬ 
tion  to  the  Empirical  Data 

The  next  and  the  last  step  in  the  non-parametric  estimation  of  the  transfer 
function  is  to  actually  obtain  a  finite  dimensional  transfer  function  that  describes 
the  input-output  behavior  of  the  system  with  a  reasonable  accuracy.  Inspecting 
the  Bode  plot  shown  in  Fig.  (6.11),  one  note  that  it  is  well  localized  in  the 
frequency  domain.  Roughly  speaking,  this  means  that  the  transfer  function  is 
more  concentrated  at  certain  frequencies. 
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e 

7 

a0 

bo 

1.0 

1.60 

2.0 

3.33 

Table  6.1:  Parameters  of  the  analyzing  wavelet  and  translation-dilation  step 
sizes 

The  shape  of  the  transfer  function  (localization  property)  together  with  the 
stability  of  the  overall  system  (due  to  internal  and  viscous  damping)  motivate  us 
to  use  the  method  of  rational  wavelets  to  approximate  the  transfer  function  [2]. 

This  method  has  shown  to  be  superior  (e.g.  in  terms  of  the  degree  of  the 
transfer  function  used  as  approximation)  to  the  conventional  methods  such  as 
the  Laguerre  model. 

There  are  several  parameters  that  one  can  change  to  get  a  good  approxi¬ 
mation  with  rational  wavelets.  One  of  the  important  factors  is  the  analyzing 
wavelet.  Unfortunately,  there  is  no  systematic  way  to  choose  such  a  wavelet  and 
usually  one  has  to  find  a  proper  analyzing  wavelet  through  experiment. 

The  translation  step  size  is  another  parameter  that  has  a  dramatic  effect  on 
the  approximation.  In  Fig.  (6.12)  the  result  of  a  high  order  approximation  is 
shown.  For  practical  purposes  we  have  only  considered  the  frequency  responce 
up  to  200 rad/ sec.  Figure  (6.13)  shows  the  results  for  a  rational  transfer  function 
of  degree  42. 

The  following  function  is  used  as  the  analyzing  wavelet 


= 


1 

(s-f  7)2  +  C 


(6.8.1) 


The  dilation  and  translation  step  sizes,  7,  and  £  are  as  in  Table  (6.1). 
Figure  6.14  depicts  the  frequency-time  localization  of  the  empirical  transfer 
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Figure  6.12:  High  order  approximation  (least  squares) 

function  obtained.  It  can  be  seen  from  this  figure  that  the  certain  coefficients 
are  considerably  bigger  (and  thus  more  important)  than  the  other  and  the  usage 
of  wavelet  systems  for  approximation  is  justified.  Also  the  poles  and  zeros  of 
the  approximated  transfer  function  are  shown  in  6.15. 


6.9  Closed  Form  of  the  Transfer  Function 

In  order  to  compare  the  empirical  transfer  function  obtained  in  the  previous 
section  with  the  analytical  one,  here  we  derive  a  closed  form  expression  for  the 
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Figure  6.13:  Approximation  with  a  transfer  function  of  degree  42  (least  squares) 
transfer  function. 

This  transfer  function  would  not  be  a  rational  one  because  it  is  derived  from 
a  partial  differential  equation  rather  than  an  ODE.  Consider  the  following  PDE 
which  corresponds  to  the  vibration  without  damping, 

PUtt{x,t)  +  Mxx(x,t)  =  0, 

M(x,t)  —  EIyxx(x,t)  +  CpV(x,t), 

2/(0,  f)  =  yx(0,t)  =  0,  (6.9.1) 

y{x,0)  =  yt(x,  0)  =  0, 

M(L,t)  =  Mx(L,t)  =  0. 
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Dilations 


Translations 


4 

Figure  6.14:  Magnitudes  of  the  coefficients  of  the  wavelet  system  vs.  dilation 
and  translation 

Take  the  Laplace  transform  (with  respect  to  t)  to  eliminate  t  and  obtain  a 
differential  equation  with  respect  to  x. 


ps2Y(x,s)  +  ^M{x,s)  =  0, 

M(x,s)  =  El£,Y(x1s)  +  CpV(x,s ), 

< 

y(M)  =  =  °> 

M(L,s)  =  Mx{L,s )  =  0. 

Define  the  state  variables  as, 


(6.9.2) 
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Figure  6.15:  Pole  and  zeroes  of  the  approximated  transfer  function 


ux(a:) 

= 

u2{x ) 

= 

u3(x) 

=  M(x,s), 

u^(x) 

= 

ui(0) 

=  M2(0)  =  0, 

U3(L) 

=  u4(L)  =  0. 

The  state  space  representation  would  be 


(6.9.3) 


U\x)  = 


0 

1 

0 

0 

0 

0 

0 

1 

El 

0 

U{x)  + 

-cp 

El 

0 

0 

0 

1 

0 

— ps2 

0 

0 

0 

0 

(6.9.4) 
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The  only  difficulty  with  the  above  ODE  is  that  ut  and  u2  are  specified  at 
one  end  while  u3  and  u4  are  given  at  the  other  end.  The  variation  of  constant 
for  the  system  (6.9.4)  is 


U(x) 


iAxU{x)  +  J*  eA^x-^BV((T,  s)d<r, 


(6.9.5) 


with  the  obvious  definitions  for  the  matrices  A  and  B  and  the  initial  condition, 


m 


0 

0 

u3(0) 

u4(0) 


(6.9.6) 


To  obtain  u3(0)  and  u4(0)  evaluate  Eq.  (6.9.5)  at  x  —  L  and  collect  the  last 
two  rows.  i.e. 


=  I<eALu{0)  + KiT  eA{L~a)Bda)F(s),  (6.9.7) 

J  X\ 

where  we  assumed  a  uniform  electric  field,  with  respect  to  x,  between  x\  and  x2 
and  zero  elsewhere  and  with  the  Laplace  transform  F(s)  (see  Fig.  (6.16).  Also 
K  is  given  by, 


0  = 


U3  (L) 
u4(L) 


K 


0  0  10 
0  0  0  1 


(6.9.8) 


Solving  for  u3(0)  and  u4(0)  we  get, 


u3{0) 

=  -{KeALKT)~lK 

'  reML-ABda 

u4(0) 

(6.9.9) 
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Figure  6.16:  The  electric  field  across  the  piezoceramic  at  one  instant  of  time 


Since  we  are  only  interested  in  deriving  the  transfer  function  and  the  output 
is  proportional  to  the  difference  in  slope  of  the  beam  (Eq.  (6.7.9)),  only  the 
second  state  variable  is  important.  Thus,  to  compute  u2(x)  substitute  for  the 
initial  condition  into  Eq.  (6.9.5)  and  pick  the  second  row. 


bout 


0  K,  0  0 


(Us2-USl), 


(6.9.10) 


where,  Si  and  S2  specify  the  position  of  the  sensor  and  Ks  is  as  in  Eq.  (6.9.5). 
Using  the  above  equations  and  with  the  use  of  Mathematica  the  following  ex¬ 
pression  has  been  obtained  for  the  transfer  function, 


Vout{s)  N(s) 
F{s)  ~  D(sy 

where,  N(s)  and  D(s)  are  defined  as, 


(6.9.11) 
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N(s)  =  {^^[2cos(32/4)  cosh(32/4)  —  2cos(2:/2)  cosh(z/2)] 

[—  cosh(3r:/4)  sin(3z/4)  -f  cosh(3z/4)  sin(5,z/4)— 
cos(3z/4)  sinh(3z/4)  +  cos(3z/4)  sinh(5;z/4)]  + 

[cosh(z/2)  sin(z/2)  —  cosh(32,/4)  sin(3z/4)— 
cos(z/2)  sinh(,z/2)  +  cos(3z/4)  sinh(3z/4)] 

[cos(5z/4)  cosh(32,/4)  —  cos(3,z/4)  cosh(5-z/4)+ 

2sin(3z/4)  sinh(3,z/4)  —  sin(5.z/4)  sinh(3.z/4)  —  sin(3z/4)  sinh(5.z/4)]}, 
D(s )  =  2  +  cos(z)  +  cosh(.z), 


•*  =  V2L(^Vs- 

To  obtain  the  natural  frequencies,  one  can  replace  s  with  jto  and  after  sim¬ 
plification  we  have, 

1+  cos  cosh  =  0.  (6.9.12) 

The  magnitude  of  the  transfer  function  versus  frequency  is  shown  in  Fig. 
6.17. 

One  can  follow  the  exact  same  procedure  to  get  an  expression  for  the  transfer 
function  without  ignoring  passive  damping  by  replacing  the  A  matrix  in  Eq. 
(6.9.4)  with, 

1  0  0 

CDl-ti  Wi  0 

0  0  1 

0  0  0 


0 

0 

0 

— ps 2  —  7  s 
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Figure  6.17:  Modulus  of  the  closed  form  transfer  function  vs.  frequency 

However  the  expression  in  that  case  would  be  very  lengthy.  Again  it  can  be 
verified  that  the  natural  frequencies  computed  in  Eq.  (6.7.11)  are  actually  roots 
of  the  denominator  of  the  transfer  function  obtained  (Eq.  (6.9.12))  and  also  the 
peaks  of  the  magnitude  plot  in  Fig.  (6.17)  occurs  almost  at  the  same  frequencies 
as  in  Fig.  (6.11)  which  shows  that  the  estimation  method  for  computing  the 
empirical  transfer  function  is  acceptable.  Of  course  these  two  graphs  are  different 
because  here  for  the  sake  of  simplicity  passive  damping  has  not  been  included. 
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Chapter  7 


Conclusion  and  Suggestions  for  Future 
Work 


In  this  work  we  consider  the  problem  of  active  vibration  damping  in  flexible 
structures  using  piezo  ceramics  and  some  control  strategies  are  tested  using  a 
simulation  program  written  for  this  purpose. 

Specifically,  the  application  of  a  Gutman-Hagander  type  controller  is  shown 
to  have  both  the  advantage  of  stability  under  the  assumption  that  there  is  a 
limit  on  the  input  signals  that  can  be  applied  to  the  system  and  a  faster  settling 
time  compared  to  the  conventional  LQR  controller  design. 

Although  the  performance  of  these  control  methods  has  been  tested  through 
simulation  programs,  the  actual  implementation  of  these  control  designs  into  a 
cantilever  beam  is  an  interesting  project  that  has  not  been  done  in  the  presented 
work. 

In  the  last  chapter,  we  developed  an  algorithm  for  nonparametric  system 
identification  and  as  a  particular  application  we  considered  the  problem  of  fitting 
a  transfer  function  with  piezo  ceramics  acting  both,  as  sensors  and  actuators. 

Through  the  usage  of  this  method,  one  can  obtain  the  transfer  function  (an 
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approximation)  from  each  input  to  each  output  and  then  apply  a  MIMO  control 
strategy  to  do  the  task  of  vibration  damping.  However,  the  controller  design 
problem  remains  as  a  topic  for  further  research  and  is  not  addressed  in  this 
thesis. 

Also,  once  the  controller  is  designed,  the  performance  of  it  may  be  tested  by 
using  either  the  simulation  program  written  for  this  purpose  or  by  the  actual 
implementation  of  it  in  an  experiment. 

Throughout  this  work,  we  did  not  consider  any  kind  of  distributed  parameter 
control  method  for  solving  the  problem  at  hand  (which  may  seem  to  be  more 
natural)  and  this  again  may  be  a  challenging  research  problem. 

For  a  general  treatment  of  the  distributed  parameter  control  one  may  see 
the  following  references  :  [16],  [7],  [12],  and,  [10]. 
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